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INTRODUCTION 


This  postdoctoral  traineeship  grant  (W81XWH-06-1-0235,  entitled  “Automated  patient 
positioning  guided  by  cone-beam  CT  (CBCT)  for  prostate  radiotherapy”)  was  awarded  to 
Dr.  Tianfang  Li,  the  original  principal  investigator  (PI)  on  January  1st,  2006.  Dr.  Tianfang 
Li  left  Stanford  University  and  accepted  a  faculty  position  at  the  University  of  Texas 
Southwestern  Medical  Center  at  Dallas  in  the  end  of  2006.  The  award  was  transferred  to 
the  current  PI,  Ming  Chao  Ph.D.  It  took  a  few  months  for  the  transition  paper  work  to 
complete.  A  no-cost-extension  was  filed  to  extend  the  period  of  the  project  for  one  year, 
namely,  through  December  31st,  2008.  This  is  the  final  report  for  the  entire  funding 
period  (January  1st,  2006  -  December  31st,  2008). 

The  goal  of  this  project  is  to  develop  a  clinically  practical  technique  for  prostate 
patient  positioning  based  on  newly  emerged  CBCT  on-board  imaging  system.  Under  the 
generous  support  from  the  U.S.  Army  Medical  Research  and  Materiel  Command 
(USAMRMC),  the  Pis  have  gained  a  tremendous  amount  of  knowledge  on  prostate 
cancer  and  prostate  cancer  management.  The  support  has  also  made  it  possible  for  the  Pis 
to  contribute  significantly  to  prostate  cancer  research.  A  number  of  conference  abstracts 
and  peer-reviewed  publications  have  been  resulted  from  the  support.  In  this  report,  the 
Pis’  research  activities  in  the  funding  years  are  summarized  and  the  accomplishments  of 
the  proposed  tasks  in  the  Statement  of  Work  are  presented. 


BODY 


I.  Introduction 

Modem  radiotherapy  equipment  is  capable  of  delivering  high  precision  conformal  dose 
distributions  to  the  target.  However,  how  to  precisely  locate  the  target  prior  to  each 
treatment  fraction,  especially  for  soft  tissue,  remains  a  problem,  due  to  possible  non-rigid 
internal  motion  relative  to  bony  structures  or  external  landmarks.  For  this  reason,  gold- 
seed  implants  have  been  often  used,  for  example  in 
prostate,  to  enhance  the  visibility  of  the  target  in  the 
portal  film  or  EPID  (electronic  portal  image  device) 
image.  Recently,  a  new  technology  of  kilo-voltage 
cone-beam  CT  (CBCT)  has  been  integrated  onboard 
with  the  linear  accelerator  treatment  machine.  Figure  1 
shows  a  CBCT  on-board  imager  on  a  Varian  Trilogy™ 
system.  Superior  to  the  common  approach  based  on  the 
two  orthogonal  images  provided  by  the  mega- voltage  Figure  1.  Varian  Trilogy  oncology 

EPID,  CBCT  can  provide  high-resolution  three-  system  with  on-board  CBCT  imager 

dimensional  (3D)  information  of  the  prostate  and  other 

pelvic  structures1.  Thus  the  prostate  target  localization  could  be  potentially  done  more 
accurately  with  the  aid  of  the  on-board  CBCT  imager,  even  without  using  any  implanted 
fiducials.  However,  in  practice,  there  is  currently  still  a  general  lack  of  efficient  method 
to  utilize  the  3D  CBCT  images  for  accurate  patient  positioning.  The  two  major  inherent 
difficulties  in  the  task  of  prostate  localization  are:  (1)  the  inter- fraction  organ  motions  are 
more  pronounced  at  pelvis  than  other  sites,  due  to  the  involuntary  physiological  motions 
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of  the  involved  organs  (such  as  the  rectum  and  bladder)2, 3  (Appendices  5  and  6);  (2)  the 
prostate  has  relatively  low  contrast  resolution  in  CT  images,  which  is  likely  to  cause 
more  errors  in  target  delineation  or  localization4  (Appendices  1  and  2  ).  Extensive  efforts 
have  been  made  in  utilizing  CBCT  for  prostate  target  positioning,  for  instance,  to  position 
the  patient  using  bony  landmarks.  Since  the  prostate  gland  can  move  independently  to  the 
bony  structures  and  its  shape  may  change  from  time  to  time,  this  commonly  used 
positioning  technique  results  in  large  uncertainties  in  prostate  targeting.  Manually 
matching  the  prostate  target  is  possible  but  it  is  very  time  consuming  and  will  introduce 
significant  inter/intra-operator  uncertainties.  The  goal  of  this  research  project  is  to 
develop  effective  volumetric  image  guidance  techniques  to  facilitate  the  patient 
positioning  for  prostate  radiation  therapy. 

II.  New  Patient  Positioning  Technique 

We  have  developed  a  novel  automated  patient  positioning  technique  for  prostate 
cancer  patients  using  deformable  image  registration.  The  significant  advantage  of  this 
technique  over  traditional  method  is  two-fold.  First,  it  accounts  for  the  complex  non-rigid 
motion  of  the  soft-tissue  organs;  Secondly,  it  considers  both  the  prostate  target  and  the 
organs  at  risk  (OARs)  simultaneously,  thus  a  better  clinical  decision  can  be  made  in  term 
of  balancing  the  displacements  of  the  target  and  OARs.  In  brief,  the  method  first  registers 
the  3D  CBCT  image  acquired  in  the  treatment  room,  to  the  3D  planning  CT  image  using 
B-Spline  deformable  model,  so  that  a  point-to-point  correspondence  for  the  patient  over 
the  time  (the  period  from  planning  to  treatment)  or  the  “deformation  field”  can  be 
established.  The  displacement  of  every  point  in  the  target  and  sensitive  structures  are 
subsequently  calculated  and  used  to  determine  the  final  couch  shifts  by  a  weighted  sum1. 

TM 

We  have  verified  our  positioning  technique  using  a  Varian  Trilogy  system  (Varian 
Medical  Systems,  Palo  Alto,  CA)  installed  in  our  clinic,  which  is  equipped  with  a  kilo- 
voltage  CBCT  imager  (see  Figure  1).  To  test  the  technique,  we  have  also  built  a  pelvis 
phantom  consisting  of  bladder,  prostate,  rectum,  and  bony  structures,  as  illustrated  in 
Figure  2.  The  phantom  was  first  scanned  with  a  clinical  CT  (General  Electric  Medical 
Systems,  Waukesha,  WI),  with  which 
a  four-field  3D  conformal  treatment 
plan  was  generated.  After  deformation, 
it  was  scanned  again  with  both  the 
regular  CT  and  a  The  planning  CT  was 
acquired  at  120  kV,  380  mA,  2.5mm 
slice  thickness  using  helical  mode, 
while  the  positioning  CT  (and  CBCT) 
were  acquired  at  80  mA  with  other 
scan  parameters  kept  the  same.  The 
position  of  the  “prostate”  target  was 
then  derived  using  these  images  by 
three  different  approaches,  including 
the  proposed  one  based  on  deformable 
registration,  one  based  on  bony 
anatomy  alignment,  and  one  based  on  minimizing  CT-number  difference  (MCD)  via  rigid 
translation  of  the  target5. 


Figure  2.  Schematic  Drawing  of  a  deformable 
pelvic  phantom. 
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It  was  found  that  both  the  proposed  and  the  MCD  methods  showed  clear  differences 
from  bony  anatomy  alignment,  reflecting  the  fact  that  the  target  had  an  internal  motion 
independent  to  the  bony  structures.  However,  significant  discrepancies  were  also 
observed  between  the  proposed  and  the  MCD  methods,  specifically,  12mm  in  x,  9mm  in 
y,  and  5mm  in  z  axis  in  this  study.  Comparing  Figures  3(a)  and  3(b),  we  found  that  the 
proposed  method  presented  much  better  accuracy  for  the  target  localization.  In  each  of 
the  two  figures,  positioning  parameters  were  also  calculated  by  using  the  regular  CT 
images  of  the  deformed  phantom  instead  of  CBCT,  in  which  consistent  results  were 
achieved  by  the  proposed  one  (see  Figure  3b),  while  the  MCD  method  generated  different 
parameters  (see  Figure  3a).  One  of  the  reasons  for  MCD  generating  varied  results  when 
CBCT  was  replaced  by  clinical  CT  for  positioning,  is  because  that  the  clinical  CT  and  the 
onboard  CBCT  have  observable  CT-number  difference  for  the  same  tissue. 


Figure  3.  Positioning  with  (a)  minimum  CT-number  difference  and  (b)  deformable  registration.  The  box  outlines 
the  volumetric  image  of  the  phantom  in  the  treatment  position.  The  color  contour  shows  the  same  target  after  CBCT- 
guided  repositioning,  and  the  gray  contour  is  the  target  after  regular  CT-guided  repositioning.  The  CT-number 
difference  between  clinical  CT  and  on-board  CBCT  resulted  in  discrepancy  in  the  positioning  parameters  by  this 
method. 


The  proposed  registration-based  CBCT  repositioning  is  quite  robust  to  noise  and 
system  errors  of  the  imaging  modalities.  With  this  technique,  the  uncertainties  of  soft- 
tissue  target  localization  could  be  greatly  reduced,  which  ensures  conformal  dose 
distribution  to  be  precisely  delivered  as  planned.  This  will  also  enable  the  implementation 
of  conformal  radiotherapy  for  prostate  tumor  with  much  smaller  margins  than  currently 
applied,  leading  to  less  complications  and  side  effects  and  improved  outcome  of 
radiotherapy.  The  technique  provides  an  outstanding  opportunity  for  prostate  dose 
escalation  study. 

III.  Registration  Technique  Refinement 
III.A  Theoretical  investigation 

One  of  the  most  important  issues  of  using  CBCT  images  is  to  accurately  register  it 
with  the  planning  CT  image,  which  is  not  only  essential  to  patient  positioning,  but  also 
critical  to  various  applications  such  as  treatment  planning  and  delivered  dose  verification 
or  adaptive  radiotherapy.  Currently,  the  CBCT  imaging  system  is  far  from  being  perfect. 
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Many  image  degrading  factors  such  as  the  scatters,  the  beam  hardening,  and  the 
approximate  reconstruction  algorithm  etc,  adversely  affect  the  final  image  quality  and 
may  subsequently  lead  to  inaccurate  registration.  Generally  speaking,  the  CBCT  3D 
image  is  an  artifacts-contaminated  image  compared  with  the  diagnostic  CT  image.  To 
improve  the  registration  accuracy,  we  have  been  investigating  a  novel  image  space-to- 
projection  space  registration  technique.  The  method  has  been  validated  with  computer 
simulation,  and  the  clinical  study  is  under  progress. 

Since  the  reconstruction  process  often  complicates  the  property  of  various  artifacts, 
we  may  avoid  the  drawback  of  using  low-quality  reconstructed  image  in  deriving  the 
motion  model  by  registering  the  3D  planning  CT  directly  to  the  CBCT  raw  data.  The 
unique  problem  with  this  method  is  that  the  two  images  being  registered  are  in  different 
spaces,  i.e.  one  in  the  image  domain,  and  the  other  in  the  projection  domain.  We  have 
developed  an  efficient  algorithm  for  such  a  registration  task  using  a  BSpline  model,  and 
one  example  of  the  advantage  of  such  method  is  shown  in  Figure  4. 


Figure  4.  Deformation  field  resulted  from  different  registration  technique,  (a)  direct  CT  to  CBCT  image 
registration,  where  view-aliasing  artifacts  causes  the  registration  error;  (b)  CT  image  to  CBCT  raw  data  registration;  (c) 
true  deformation  field. 


III.B  Refinement  of  image  registration  technique 

The  novel  automated  patient  positioning  technique  was  developed  for  prostate  cancer 
patients  using  deformable  image  registration.  Compared  to  the  conventional  approach, 
the  new  technique  accounts  for  the  complex  non-rigid  motion  of  the  soft-tissue  organs 
and  balances  the  prostate  target  and  the  organs  at  risk  simultaneously.  The  method  was 
verified  using  a  pelvis  phantom  consisting  of  bladder,  prostate,  rectum  and  bony 
structures.  We  have  further  refined  the  registration  technique  and  significantly  improved 
the  computational  efficiency  of  the  technique,  thus  making  it  possible  for  future  clinical 
application.  Briefly,  our  efforts  were  focused  on  improving  the  registration  technique  by 
using  different  similarity  measures  such  as  normalized  cross  correlation  and  mutual 
information  (Appendices  1,  3,  and  4).  We  found  that  by  using  the  Mattes  Mutual 
Information  (MMI)  metric,  not  only  significant  amount  of  computing  time  was  saved,  but 
also  precision  of  the  registration  was  achieved2.  The  central  concept  of  mutual 
information  (MI)  is  the  calculation  of  entropy.  For  an  image  A,  the  entropy  is  defined  as 

H (A)  = -^ p A(a)\og p A{a)da  ,  where  pA (a)  (also  termed  as  marginal  probability 
density  function  (PDF))  is  the  probability  distribution  of  grey  values  (image  intensities) 
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which  is  estimated  by  counting  the  number  of  times  each  grey  value  occurs  in  the  image 
and  dividing  those  numbers  by  the  total  number  of  occurrences.  Given  two  images,  A  and 

B,  their  joint  entropy  is  H(A,B)  =  -^  pAB(a,b)\ogpAB(a,b)dadb,  wher epAB(a,b)  is 

the  joint  PDF  defined  by  a  ratio  between  the  number  of  grey  values  in  the  joint  histogram 
(feature  space)  of  two  images  and  the  total  entries.  The  mutual  information  is  generally 
expressed  as  MI  {A,  B )  =  H  (A)  +  H(B )  -  H(A,  B) .  MI  measures  the  level  of  information 
that  a  random  variable  (e.g.,  Ia(x))  can  predict  about  another  random  variable  (e.g.  Ib(x)). 
Different  from  the  conventional  MI,  where  two  separate  intensity  samples  are  drawn 
from  the  image,  the  Mattes  implementation,  MMI,  uses  only  one  set  of  intensity  to 
evaluate  both  the  marginal  and  joint  PDFs  at  discrete  positions  or  bins  that  uniformly 
spread  within  the  dynamic  range  of  the  images6.  Entropy  values  were  computed  by 
summing  over  all  the  bins.  The  introduction  of  the  mutual  information  metric  has  showed 
significant  improvement  on  the  registration  precision  and  helped  reducing  the  impact  of 
scatter  and  other  factors  inherent  in  CBCT  imaging. 

IV.  Novel  strategies  of  enhancing  the  quality  of  CBCT  images 

Onboard  CBCT  imaging  plays  an  essential  role  in  the  success  of  prostate  IMRT  ' . 
Unfortunately,  the  image  quality  of  currently  available  CBCT  is  far  from  satisfactory  due  to  a 
number  of  issues  specific  to  the  CBCT  scan  and  image  reconstruction.  Factors  that  adversely 
affect  the  CBCT  performance  include  scatter  photons,  motion  artifacts,  and  excessive  radiation 
dose.  We  were  among  the  first  in  developing  method  to  utilize  the  prior  knowledge  of  the 
system  to  improve  the  resultant  CBCT  images.  Generally,  the  patient  has  already  had  a 
planning  CT  (and  possibly  one  or  more  CBCTs)  before  an  on-treatment  CBCT  is  acquired. 
These  imaging  data  can  be  utilized  as  a  priori  knowledge  to  enhance  the  image  quality  and  to 
reduce  the  radiation  dose  in  routine  CBCT  imaging  process.  This  has  been  successfully 
demonstrated  in  our  previous  studies  on  CT  and  PET  imaging10.  For  CBCT  imaging,  due  to  the 
significantly  reduced  number  of  projections  per  reconstruction,  the  quality  of  the  phase- 
resolved  CBCT  images  is  greatly  degraded  by  the  view-aliasing  artifacts.  Most  commonly, 
acquisitions  using  multiple  gantry  rotations  or  slow  gantry  rotation  can  increase  the  number  of 
projections  and  subsequently  improve  the  image  quality.  However,  the  price  to  pay  here  is  the 
significant  increase  of  scanning  time.  Our  CBCT  image  enhancement  method  effectively 
circumvent  this  problem  through  effective  utilization  of  the  information  contained  in  other 
phases  and/or  even  previously  obtained  planning  CT  (fan  beam  CT).  By  registering  the  prior 
data  to  the  CBCT  (in  either  image  space  or  projection  space),  we  were  able  to  enhance  the 
contrast-to-noise  ratio  (CNR)  of  the  CBCT  images  by  5~10  fold,  without  noticeable 
compromise  in  the  spatial  resolution  of  the  resultant  images.  The  developed  technique 
significantly  reduces  the  view-aliasing  artifacts  in  CBCT  images  and  leads  to  high  quality 
images  for  therapeutic  guidance.  This  technique  also  dramatically  decreases  the  radiation  dose 
of  CBCT  acquisition.  The  image  quality  enhancement  makes  it  possible  to  directly  use  CBCT 
for  accurate  dose  calculation  and  IMRT  planning. 

In  summary,  along  the  direction  of  Task  4  in  the  Statement  of  Work,  we  have  done 
substantial  work  on  improving  the  performance  of  the  CBCT  system.  We  developed  methods 
for  enhancing  the  performance  of  CBCT  and  reducing  radiation  dose  incurred  during  frequent 
CBCT  scans4,  n.  These  accomplishments  are  attached  in  the  Appendices  (items  1  and  2),  which 
conclude  the  task  4  in  the  Statement  of  Work. 


8 


V.  Exploration  of  Low-Dose  CBCT  Imaging 

Due  to  the  repeated  daily  use  of  CBCT  for  a  patient,  the  radiation  exposure  of  the 
imaging  system  has  been  of  great  concern.  For  patient  positioning  purpose,  the  strategies 
to  reduce  the  CBCT  radiation  exposure  that  we  have  developed  during  this  training 
period  are  using  the  lowest  possible  x-ray  tube  current  and  low  number  of  projections. 
The  issues  associated  with  low  x-ray  tube  current  is  mainly  the  dramatically  increased 
statistical  noise,  and  with  low  number  of  projection,  the  undersampling  (view-aliasing) 
artifacts  become  prominent  in  CBCT  images.  Apparently,  simple  reduction  of  the 
radiation  dose  results  in  bad  image  quality.  In  order  to  maintain  an  acceptable  image 
quality  or  a  minimal  patient  positioning  error,  we  have  proposed  general  CBCT 
enhancing  techniques,  which  are  based  on  the  deformable  registration  of  all  available 
information  (CT  or  CBCT  images)  of  the  same  patient  at  different  time.  The  extra 
information  borrowed  from  other  images  helps  to  improve  the  statistics,  and  the  data 
sampling  rate  as  well.  As  one  of  the  examples  we  have  studied,  in  Figure  5  it  is  shown  the 
improvement  of  our  enhancing  technique  with  low-dose  CBCT  acquisition  strategy. 


Figure5.  Left:  CBCT  image  with  low  number  of  projections  showing  under-sampling  streak  artifacts;  Right:  the 
improved  CBCT  image  with  our  enhancing  technique  based  on  deformable  registration. 


17.  Contour  propagation  from  planning  CT  to  CBCT 

Large  interfraction  organ  motion  exists  in  radiation  therapy  of  pelvic  diseases, 
such  as  prostate  and  cervical  cancers.  Adaptive  therapy  provides  a  viable  option  to  ensure 
an  adequate  dose  coverage  of  the  tumor  target  while  sparing  the  sensitive  structures12, 13. 
Toward  the  general  goal  of  more  accurate  patient  positioning  and  clinically 
implementation  of  adaptive  therapy  of  prostate  cancer,  we  developed  an  effective 
regional  contour  mapping  technique  to  automatically  propagate  rectum  and  prostate 
contours  from  planning  CT  to  CBCT2, 14-16  (Appendices  3,  4,  5,  and  6).  This  effort  is 
based  on  the  precise  deformable  registration  between  planning  CT  and  CBCT  images. 
Different  from  disease  sites,  such  as  the  lungs  and  liver,  the  contour  mapping  here  in 
prostate  radiation  therapy  is  complicated  by  two  factors:  (i)  the  physical  one-to-one 
correspondence  may  not  exist  due  to  the  insertion  or  removal  of  some  image  contents 
within  the  rectum;  and  (ii)  reduced  signal  to  contrast  ration  of  the  CBCT  images  due  to 
increased  scatter  in  CBCT  scanning.  A  solution  customized  to  rectum  mapping  is 
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proposed  to  overcome  the  above  two  limitations  and  allow  us  to  take  advantage  of  the 
regional  calculation  algorithm  yet  avoiding  the  nuisance  of  rectum/bladder  filling.  The 
implementation  of  the  contour  propagation  was  achieved  with  the  narrow  shell  technique. 
The  shell  was  constructed  outside  the  manually  delineated  contours  on  planning  CT 
including  the  region  of  1  ~  2  cm.  The  central  idea  here  is  to  exclude  the  volume  inside  the 
rectum  because  its  content  may  change  from  day  to  day.  Prostate  contour  can  be  mapped 
similarly,  but  the  shell  spans  the  regions  both  inside  and  outside  the  manually  segmented 
prostate.  Figure  6  illustrates  the  planning  CT  and  the  template  contours  for  the  ROIs  such 
as  prostate  and  rectum.  Based  on  the  narrow  shell  technique,  various  deformable  models 
such  as  B-Spline  and  Thin  Plate  Spline  models  were  employed  to  achieve  the  goal  of 
contour  propagation  ’  .  A  feature  based  model  using  the  Scale  Invariant  Feature 

Transformation  (SIFT)  was  also  explored  for  higher  registration  accuracy16,  19,  20 
(Appendix  5).  The  mapped  rectal  contour  is  illustrated  in  Figure  7.  This  task  was 
completed  during  the  second  funding  period  -  one  manuscript  reporting  the  use  of  narrow 
shell  based  contour  mapping  strategy  has  been  accepted  by  Physics  in  Medicine  and 
Biology  for  publication.  In  addition,  a  paper  reporting  the  use  of  SIFT  method  to  auto- 
identify  the  tissue  features  for  robust  deformable  registration  and  contour  mapping  in  the 
pelvic  region  has  been  conditionally  accepted  by  Medical  Physics  for  publication 
(Appendices  5  and  6). 


Figure  6.  Narrow  shell  representation 
of  regions  of  interest. 


Figure  7.  Mapped  rectal  contours  (in 
red).  Manual  contour  is  overlaid  (blue) 


KEY  RESEARCH  ACCOMPLISHMENTS 

•  Developed  a  new  image-to-projection  deformable  registration  algorithm  for 
CBCT/CT  image  matching,  which  is  more  accurate  than  conventional  image-to- 
image  registration  method  and  more  robust  to  various  image  artifacts  such  as  view¬ 
aliasing  or  truncations  (Task  1  of  SOW). 

•  Developed  and  clinically  implemented  a  novel  automatic  patient  positioning  strategy, 
which  accounts  for  non-rigid  organ  motion,  and  balance  the  displacements  of  various 
organ  structures  (Task  2  of  SOW). 
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•  Improved  the  CBCT  imaging  quality  and  therefore  made  it  possible  to  use  low 
radiation  dose  CBCT  for  daily  use,  significantly  reduced  the  risk  of  radiation-induced 
secondary  cancer  for  patients  (Task  4  of  SOW). 

•  Established  a  novel  technique  to  enhance  on-board  CBCT  and  to  effectively  reduce 
the  radiation  dose  (Task  4  of  SOW) 

•  Improved  the  registration  technique  by  investigating  various  similarity  measures 
including  both  the  normalized  cross  correlation  metric  and  the  Mattes  mutual 
information  metric  (Task  2  of  SOW) 

•  Developed  a  novel  technique  of  narrow  shell  representation  of  the  regions  of  interest 
better  target  localization  and  critical  structure  modeling  to  account  for  the  large  inter¬ 
fraction  (Task  3  of  SOW) 

•  Established  the  framework  of  contour  propagation  from  planning  CT  to  CBCT  based 
on  the  narrow  shell  technique  for  the  adaptive  prostate  radiation  therapy  (Task  3  of 
SOW) 
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publication  materials  are  enclosed  with  this  report. 
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•  T.  Li,  A.  Koong,  L.  Xing,  “Enhanced  4D  cone-beam  CT  with  inter-phase  motion  model”, 
Medical  Physics,  34  3688-3695,  2007 

•  T.  Li,  L.  Xing,  “Optimizing  4D  cone-beam  CT  acquisition  protocol  for  external  beam 
radiotherapy”,  International  Journal  of  Radiation  Oncology,  Biology,  Physics,  67  1211- 
1219,2007 

•  M.  Chao,  E.  Schreibmann,  T.  Li,  N.  Wink,  L.  Xing,  “Automated  contour  mapping  using 
sparse  volume  sampling  for  4D  radiation  therapy”,  Medical  Physics,  34  4023-4029,  2007 

•  N.  Wink,  M.  Chao,  J.  Antony,  L.  Xing,  “Individualize  treatment  margin  in  respiration¬ 
gated  radiation  therapy”,  Physics  in  Medicine  and  Biology,  53  165-175,  2008 

•  M.  Chao,  T.  Li,  E.  Schreibmann,  A.  Koong,  L.  Xing,  “Automated  contour  mapping  with 
a  regional  deformable  model”,  International  Journal  of  Radiation  Oncology,  Biology, 
Physics,  70  599-608,  2008 

•  Y.  Xie,  M.  Chao,  P.  Lee,  L.  Xing,  “Feature-based  rectal  contour  propagation  from 
planning  CT  to  cone  beam  CT  for  adaptive  radiotherapy”,  Medical  Physics,  35  4474-4487, 
2008 

•  M.  Chao,  Y.  Xie,  L.  Xing,  “Auto-propagation  of  contours  for  adaptive  prostate  radiation 
therapy”,  Physics  in  Medicine  and  Biology,  53  4533-4542,  2008 

Published  Abstracts 

The  Pi’s  group  has  also  been  active  in  disseminating  the  research  results.  The  following  are  some 

of  the  presentations  given  in  various  national  and  international  conferences. 

•  T  Li,  E.  Schreibmann,  A.  Koong,  Q.  Xu,  R.  Hamilton  and  L.  Xing,  “Verification  of  gated 
radiation  therapy  using  pre-treatment  4D  CBCT,”  International  Journal  of  Radiation 
Oncology  Biology  Physics,  66(3),  S604,  2006. 
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•  W.  Mao,  T.  Li,  P.  Munro,  M.  Chao,  and  L.  Xing,  “Individualizing  4D  CBCT  acquisition 
protocol  for  external  beam  radiotherapy,”  International  Journal  of  Radiation  Oncology 
Biology  Physics,  66(3),  S146-S147,  2006. 

•  L.  Xing,  A.  de  la  Zerd,  M.  Cao,  T.  Li,  B.  Armbrush,  Y.  Yang,  P.  Lee,  T.  Pawlicki,  S. 
Hancock  and  C.  King,  “On-Board  Volumetric  CT-based  Adaptive  IMRT  For  Improved 
Prostate  Cancer  Treatment,”  International  Journal  of  Radiation  Oncology  Biology 
Physics,  66(3),  S624-S625 , 2006. 

•  M.  Chao,  E.  Schreibmann,  T.  Li,  A.  Koong,  K.A.  Goodman  and  L.  Xing,  “Automatic 
Contouring  in  4D  Radiation  Therapy,”  International  Journal  of  Radiation  Oncology 
Biology  Physics,  66(3),  S649,  2006. 

•  E.  Elder,  E.  Schreibmann,  T.  Li,  T.  Fox,  L.  Xing,  J.  Crocker  and  J.  Landry,  “Registration 
of  4D  CBCT  and  4D  CT  for  Extracranial  Stereotactic  Treatments,”  International  Journal 
of  Radiation  Oncology  Biology  Physics,  66(3),  S651,  2006. 

•  P.  Peng,  M.  Chao,  Q.  Le,  T.  Li,  A.  Hsu,  T.A.  Pawlicki  and  L.  Xing,  “Auto  Contour 
Mapping  in  CBCT  for  Adaptive  Therapy  Treatment  Planning,”  International  Journal  of 
Radiation  Oncology  Biology  Physics,  66(3),  S651-S652,  2006. 

•  T.  Li,  L  Xing,  P.  Munro  et  al,  “4D  Cone-Beam  CT  (CBCT)  Using  An  On-Board  Imager,” 
Med.  Phys.  33,  2234,  2006. 

•  M.  Chao,  T.  Li,  and  L.  Xing,  “Enhanced  4D  CBCT  Imaging  for  Slow-Rotating  On-Board 
Imager  ,”  Med.  Phys.  33,  2269,  2006. 

•  M.  Chao,  E.  Schreibmann,  T.  Li,  and  L.  Xing,  “Knowledge-Based  Auto-Contouring  in 
4D  Radiation  Therapy,”  Med.  Phys.  33,  2171,  2006. 

•  M.  Chao,  E.  Schreibmann,  T.  Li,  and  L.  Xing,  “Automatic  Contour  Mapping  in  4D 
Radiation  Therapy”,  Proceedings  of  the  XVth  International  Confenrence  on  the  Use  of 
Computers  in  Radiation  Therapy  (ICCR),  June  4-6,  2007,  Toronto,  Canada 

•  C.  Wang,  M.  Chao,  L.  Xing,  “Electron  density  mapping  for  MRI-based  treatment 
planning”,  Medical  Physics  34,  2592,  2007 

•  T  Li,  L  Papiez,  R  Timmerman,  H  Choy,  A  Koong,  L  Xing,  “High-Quality  Four- 
Dimensional  CBCT  Reconstruction  with  Virtual  Projections”,  Medical  Physics  34,  2638, 
2007 

•  M.  Chao,  E.  Schreibmann,  T.  Li,  L.  Xing,  “Automated  propagation  of  region-of-interest 
contours  between  4DCT  Images  using  a  regional  deformable  model”,  Medical  Physics 
34,2515,2007 

•  N.  Wink,  M.  Chao,  L.  Xing,  “Individualized  gating  windows  based  on  four-dimensional 
CT  information  for  respiration  gated  radiotherapy”,  Medical  Physics  34,  2384,  2007 

•  M.  Chao,  L.  Xing,  “Auto-recontouring  of  CBCT  images  for  adaptive  therapy”,  Medical 
Physics  34,  2352,  2007 

•  A.  de  la  Zerd,  M.  Chao,  B.  Armbrush,  Y.  Yang,  S.  Hancock,  C.  King,  T.  Li,  L.  Xing, 
“Adaptive  IMRT  for  improved  prostate  cancer  treatment”,  Proceedings  of  the  U.S. 
Department  of  Defense  (DOD)  Prostate  Cancer  Research  Program  (PCRP)  Innovative 
Minds  in  Prostate  Cancer  Today  (IMPaCT)  meeting,  September  5-8,  2007,  Atlanta, 
Georgia 

•  M.  Chao,  Y.  Xie,  Q.  Le,  L.  Xing,  “Modeling  the  volumetric  change  of  head  and  neck 
tumor  in  response  to  radiation  therapy”,  International  Journal  of  Radiation  Oncology  * 
Biology  *  Physics,  69(3),  S741,  2007 

•  T.  La,_M.  Chao,  L.  Xing,  Q.  Le,  “Evaluation  of  intrafraction  motion  in  head  and  neck 
cancer  during  radiotherapy”,  International  Journal  of  Radiation  Oncology  *  Biology  * 
Physics ,  69(3),  S68 1  -  S682,  2007 
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•  L.  Xing,  L.  Lee,  M.  Chao,  P.  Keall,  T.  La  and  Q.  Le,  “Clinical  implementation  of  on¬ 
board  CBCT-based  adaptive  IMRT  for  head  and  neck  cancer”,  International  Journal  of 
Radiation  Oncology  *  Biology  *  Physics,  69(3),  S446,  2007 

•  J.  Wang,  M.  Chao,  L.  Xing,  “Toward  clinical  implementation  of  adaptive  treatment 
planning:  auto-propagation  of  contours  from  planning  CT  to  cone  beam  CT”, 
International  Journal  of  Radiation  Oncology  *  Biology  *  Physics,  69(3),  S43,  2007 

•  L.  Chen,  M.  Chao,  Y.  Xie,  L.  Xing,  “Tumor  shrinkage  modeling:  towards  adaptive 
radiation  therapy  for  head  and  neck  cancer”,  International  Journal  of  Radiation  Oncology 

*  Biology  *  Physics,  72(1),  S598,  2008 

•  Y.  Xie,  M.  Chao,  L.  Xing,  “Modeling  the  deformation  and  sliding  motion  of  lungs”, 
International  Journal  of  Radiation  Oncology  *  Biology  *  Physics,  72(1),  S620,  2008 

•  M.  Chao.  L.  Lee,  L.  Xing,  “Functional  lung  sparing  using  4D-CT  and  intensity- 
modulated  radiation  therapy”,  International  Journal  of  Radiation  Oncology  *  Biology  * 
Physics,  72(1),  S621,  2008 

•  J.  Antony,  G.  Luxton,  L.  Lee,  M.  Chao,  D.  Carlson,  L.  Xing,  “Biological  modeling 
indices  for  4D  radiation  therapy”,  International  Journal  of  Radiation  Oncology  *  Biology 

*  Physics,  72(1),  S629,  2008 

•  M.  Chao,  Y.  Xie,  L.  Xing,  “Image-based  modeling  of  tumor  shrinkage  or  growth: 
towards  adaptive  radiation  therapy  of  head-n-neck  cancer”,  Medical  Physics  35,  2924 
(2008) 

•  M.  Chao.  L.  Xing,  “Functional  lung  preservation  using  four-dimensional  computed 
tomography  and  intensity-modulated  radiation  therapy  for  lung  cancer”,  Medical  Physics 
35,2701,2008 

•  Y.  Xie,  M.  Chao,  L.  Xing,  “4D  CT  image-based  modeling  the  deformation  and  sliding 
motion  of  lungs”,  Medical  Physics  35,  2914,  2008 


CONCLUSIONS 

Despite  the  well-appreciated  fact  that  onboard  CBCT  system  can  provide  3D 
anatomical  information  of  a  patient  in  the  treatment  position,  its  clinical  applications  have 
been  hindered  by  the  deficiencies  in  image  quantitative  accuracy  and  by  the  lack  of  a 
comprehensive  image  analysis  procedure.  This  is  evidenced  by  that  very  few  institutions 
are  using  CBCT  routinely  for  prostate  cancer  treatment.  A  challenge  here  is  how  to 
achieve  automatic  and  accurate  patient  positioning  with  consideration  of  the  non-rigid 
motions  of  the  prostate  and  surrounding  structures. 

During  the  funding  years,  the  Pis  have  contributed  greatly  to  development  of  novel 
prostate  cancer  patient  positioning  technique.  They  have  also  gained  improved 
knowledge  on  prostate  cancer  and  prostate  radiation  therapy  treatment.  Various  factors 
affecting  the  final  positioning  decisions  are  addressed.  A  few  important  milestones  have 
been  achieved  toward  the  goal  of  the  project.  These  include:  (1)  developed  and  clinically 
implemented  an  automated  patient  positioning  strategy  and  tested  with  phantom 
experiments;  (2)  developed  image-to-projection  deformable  registration  algorithm  to 
improve  the  CBCT  image  quality  by  reducing  the  view-aliasing  artifacts;  (4)  established 
a  novel  technique  to  enhance  on-board  CBCT  and  to  effectively  reduce  the  radiation  dose 
incurred  in  the  scanning  process;  (5)  developed  a  robust  registration  model  with 
improved  metric  function;  (6)  developed  an  innovative  narrow  shell  technique  for  better 
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target  localization  and  model  the  critical  structure;  and  (7)  Established  the  method  for 
auto-propagation  of  contours  from  planning  CT  to  CBCT  based  on  feature  based  spline 
deformable  registration.  The  data  and  experiences  accumulated  in  this  project  are  well 
documented,  and  serve  as  a  good  reference  and  guidance  for  effective  use  of  volumetric 
CBCT  imaging  technology.  The  research  improves  volumetric  image  guided  prostate 
radiation  therapy  and  enhances  the  outlook  of  the  onboard  CBCT  in  clinical  practice  as  a 
whole. 

The  DOD  Postdoctoral  Training  Award  not  only  helped  the  Pis  focus  their  efforts  in 
the  field  of  radiation  therapy  of  prostate  cancer,  but  it  also  provided  them  opportunities  to 
obtain  better  understanding  and  gain  significant  experiences  in  prostate  cancer  research 
and  substantially  promoted  their  career  development  in  the  area.  Under  the  outstanding 
mentorship  of  Professor  Lei  Xing,  both  Pis  were  successful  in  their  careers.  The  initial 
PI,  Dr.  Tianfang  Li  was  offered  a  tenure  track  assistant  professor  position  at  the 
University  of  Texas  Southwestern  Medical  Center  at  Dallas  in  late  2006.  The  successor 
of  the  Award,  Dr.  Ming  Chao,  was  appointed  assistant  professorship  at  the  University  of 
Arkansas  for  Medical  Sciences  in  2008.  The  prostate  research  program  was  very 
successful  in  training  junior  researchers  and  helping  them  achieve  their  career  goal  in 
prostate  research. 
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Four-dimensional  (4D)  cone-beam  CT  (CBCT)  is  commonly  obtained  by  respiratory  phase  binning 
of  the  projections,  followed  by  independent  reconstructions  of  the  rebinned  data  in  each  phase  bin. 
Due  to  the  significantly  reduced  number  of  projections  per  reconstruction,  the  quality  of  the 
4DCBCT  images  is  often  degraded  by  view-aliasing  artifacts  easily  seen  in  the  axial  view.  Acqui¬ 
sitions  using  multiple  gantry  rotations  or  slow  gantry  rotation  can  increase  the  number  of  projec¬ 
tions  and  substantially  improve  the  4D  images.  However,  the  extra  cost  of  the  scan  time  may  set 
fundamental  limits  to  their  applications  in  clinics.  Improving  the  trade-off  between  image  quality 
and  scan  time  is  the  key  to  making  4D  onboard  imaging  practical  and  more  useful.  In  this  article, 
we  present  a  novel  technique  toward  high-quality  4DCBCT  imaging  without  prolonging  the  acqui¬ 
sition  time,  referred  to  as  the  “enhanced  4DCBCT”.  The  method  correlates  the  data  in  different 
phase  bins  and  integrates  the  internal  motion  into  the  4DCBCT  image  formulation.  Several  strate¬ 
gies  of  the  motion  derivation  are  discussed,  and  the  resultant  images  are  assessed  with  numerical 
simulations  as  well  as  a  clinical  case.  ©  2007  American  Association  of  Physicists  in  Medicine. 
[DOI:  10.1118/1.2767144] 

Key  words:  cone-beam,  4D  CT,  on-board  imager,  IGRT,  organ  motion 


I.  INTRODUCTION 


Medical  linear  accelerators  equipped  with  cone-beam  CT 
(CBCT)  imaging  system,  using  either  a  kV  or  MV  x-ray 
source,  have  recently  become  available  in  clinics.  The  volu¬ 
metric  image  provided  by  CBCT  opens  new  avenues  for  pa¬ 
tient  setup,  dose  verification,  and  online  treatment 
planning.1-12  These  applications  rely  highly  upon  the  fidelity 
and  quality  of  the  CBCT  images.  When  imaging  in  the  re¬ 
gion  of  thorax  and  upper  abdomen,  however,  respiration  in¬ 
duced  artifacts,  such  as  blurring,  doubling,  streaking,  and 
distortion  are  observed,  which  heavily  degrade  the  image 
quality  and  affect  the  target  localization  ability  and  the  accu¬ 
racy  of  the  dose  calculation.13,14 

A  recently  proposed  method  to  account  for  respiratory 
motion  during  CBCT  imaging  is  called  “respiration  corre¬ 
lated  CBCT”  or  four-dimensional  (4D)  CBCT.15-18  In  this 
method,  the  CBCT  projections  are  divided  into  several 
groups  according  to  their  respiratory  phases.  Each  phase 
group  is  then  reconstructed  independently  to  yield  a  volu¬ 
metric  image  corresponding  to  that  specific  phase.  Since  the 
selected  projections  in  each  group  have  almost  the  same 
phase,  the  method  greatly  reduces  the  motion-induced  incon¬ 
sistency  among  the  data,  leading  to  improved  image  recon¬ 
struction.  However,  the  number  of  projections  available  for 
each  image  reconstruction  is  substantially  decreased  com¬ 
pared  to  a  conventional  3D  CBCT  due  to  the  data  dividing. 
This  may  lead  to  a  serious  undersampling  problem  and  result 
in  strong  view-aliasing  artifacts  in  the  phase-resolved  im¬ 


ages.  The  artifacts  are  generated  during  the  backprojection 
step  in  the  CBCT  image  reconstruction  process  and  are  seen 
mainly  in  the  2D  axial  planes. 

To  eliminate  the  view-aliasing  artifacts  in  4DCBCT,  we 
have  recently  investigated  acquisition  techniques  of  “mul¬ 
tiple  gantry  rotation”  and  “slow  gantry  rotation”  with  low 
x-ray  tube  current.16,17  These  methods  increase  the  number 
of  projections  of  each  phase  and  significantly  improve  the 
4DCBCT  image  quality.  However  a  practical  issue  is  the 
prolonged  acquisition  time,  which  may  limit  their  clinical 
applications. 

The  purpose  of  this  work  is  to  investigate  a  novel  ap¬ 
proach  to  reduce  the  acquisition  time  while  maintaining  the 
4DCBCT  image  quality,  or  from  another  point  of  view,  to 
enhance  the  current  4DCBCT  image  quality  while  using 
short  acquisition  time.  The  approach  developed  in  this  article 
is  essentially  to  reduce  the  view-aliasing  artifacts  resulted 
from  insufficient  sampling.  To  increase  the  sampling  rate  in 
each  individual  phase,  the  projections  of  all  other  phases  can 
be  incorporated  into  the  reconstruction  process.  To  do  so,  a 
motion  model  linking  the  data  of  different  phases  is  required, 
which  can  be  derived  from  deformable  registrations  of  the 
4DCBCT  phases.  Two  other  registration  strategies  of  intro¬ 
ducing  an  additional  diagnostic  CT  image  and  mapping  be¬ 
tween  image  space  and  projection  space  are  also  investi¬ 
gated.  The  proposed  approach  is  evaluated  with  numerical 
simulations.  Since  the  view-aliasing  artifacts  are  prominent 
only  in  the  axial  view  and  for  simplicity,  we  start  with  fan- 
beam  (2D)  geometry  to  study  the  performance  of  the  ap¬ 
proach,  a  clinical  case  of  a  lung  cancer  patient  is  then  pre¬ 
sented. 
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Fig.  1.  Flow  chart  of  the  proposed  reconstruction  process  for  enhanced 
4DCBCT  imaging. 


II.  METHODS  AND  MATERIALS 

II.A.  Modified  Feldkamp  reconstruction  for  reduced 
view-aliasing  artifacts 

As  described  earlier,  to  form  the  4DCBCT  images,  the 
acquired  projection  data  are  retrospectively  sorted  into  sev¬ 
eral  phase  bins  before  reconstructions.  When  the  acquisition 
time  is  limited,  the  reduced  number  of  projections  per  recon¬ 
struction  will  result  in  undersampling  and  degrade  the  image 
with  view-aliasing  artifacts.  In  order  to  improve  the  sampling 
rate  in  each  phase,  we  need  to  borrow  projection  information 
from  other  phases.  Obviously,  direct  inclusion  of  projections 
of  other  phases  will  cause  the  inconsistency  in  the  data  due 
to  patient  motion.  However,  as  we  have  shown  previously,  if 
the  motion  can  be  accurately  modeled,  a  motion-corrected 
CBCT  reconstruction  can  be  obtained.19  This  reconstruction 
method  is  a  modified  Feldkamp  algorithm,  which  performs 
the  backprojections  of  the  filtered  data  along  deformed  paths 
according  to  the  corresponding  motion  model.  Equivalently 
and  computationally  more  efficiently,  one  can  reconstruct 
each  individual  phase  first,  then  superimpose  all  the  phases 
after  deforming  each  of  them  into  the  same  phase  based  on 
the  given  motion  model.  This  method  can  be  described  in  the 
following  how  chart  (Fig.  1). 

II. B.  Motion  model  estimation 

To  derive  the  motion  model  (i.e.,  the  points  of  correspon¬ 
dence  between  any  two  phases)  to  be  used  in  the  enhanced 
4DCBCT  reconstruction,  we  have  investigated  the  following 
three  registration  strategies  with  and  without  the  aid  of  an 
additional  3D  CT  image. 

ll.B.1.  Register  artifacts-contaminated  4DCBCT 
images 

The  most  straightforward  way  to  find  the  motion  model  is 
to  register  the  regular  4DCBCT  phases  using  a  deformable 


model.  A  unique  problem  here  is  that  the  images  to  be  reg¬ 
istered  contain  serious  steak  artifacts,  which  may  adversely 
affect  the  accuracy  of  the  derived  motion  model,  and  the 
error  can  propagate  into  the  final  4DCBCT  images.  Never¬ 
theless,  the  process  may  prove  to  be  favorable  because  of 
greatly  reduced  view-aliasing  artifacts.  The  performance  of 
this  method  is  shown  later  with  a  simulation  study  where  the 
registration  was  done  based  on  a  free-form  Spline  (B  Spline) 
model.  The  simplicity  and  yet  accuracy  of  the  B Spline 
method  make  it  a  preferred  tool  for  many  clinical 
applications.  More  details  about  this  registration  algo¬ 
rithm  can  be  found  in  Ref.  22. 


//.B.2.  Register  planning  CT  to  4DCBCT  phased 
image 

Note  that  in  radiation  therapy,  it  is  often  the  case  that  a 
patient  has  already  had  a  3D  or  4D  CT  scan  for  the  treatment 
planning  before  any  radiation  delivery.  In  such  case,  an 
artifact-free  volumetric  image  can  be  assumed  available  prior 
to  the  4DCBCT  scan  on  the  treatment  day.  Considering  res¬ 
piration  motion,  this  CT  image  should  be  one  phase  of  4D 
CT  images  or  a  breath-hold  CT  scan.  The  idea  is  to  register 
each  phase  of  the  4DCBCT  images  to  this  artifact- free  3D 
image  to  derive  corresponding  deformation  with  respect  to 
the  specific  phase  of  the  CT  image.  The  deformation  of  the 
4DCBCT  phases  with  respect  to  each  other  can  be  subse¬ 
quently  obtained  by  calculating  the  relative  movements. 
Since  only  one  of  the  two  images  being  registered  has  arti¬ 
facts,  it  is  expected  that  the  accuracy  will  be  improved  over 
the  first  procedure.  Again,  the  performance  of  this  method  is 
evaluated  with  the  simulation  study  later  in  this  article. 


II.B.3.  Registration  from  image  space  to  projection 
space 

Although  the  reconstructed  4DCBCT  phase  images  con¬ 
tain  view-aliasing  artifacts,  one  should  realize  that  the  arti¬ 
facts  are  not  originally  present  in  the  raw  data  but  are  gen¬ 
erated  later  during  the  reconstruction  process.  Therefore,  we 
may  avoid  the  drawback  of  using  a  low-quality  reconstructed 
image  in  deriving  the  motion  model  by  registering  the  3D 
planning  CT  directly  to  the  CBCT  raw  data  (after  phase  sort¬ 
ing).  Note  that  the  two  images  being  registered  are  in  differ¬ 
ent  spaces;  in  the  following,  we  present  a  simple  algorithm 
for  such  a  registration  task  using  a  B Spline  model. 

Similar  to  the  conventional  image  registration  using  de¬ 
formable  models,  the  motion  of  the  object  is  described  by 
g(x\t)=f(x+u{x\t )),  where  /( x)  represents  the  3D  image 
obtained  from  the  planning  CT  scan,  g{x\t)  is  the  same  ob¬ 
ject  at  a  particular  phase  t  during  the  4DCBCT  scan,  and 
u(x;t )  defines  the  corresponding  deformation  for  point  x  at 
that  phase.  With  the  B  Spline  model,  the  deformation  is  de¬ 
fined  on  a  sparse,  regular  grid  of  control  points  \n  placed 
over  the  CT  image.  When  a  point  is  not  on  the  grid,  the 
associated  deformation  is  calculated  by  the  B Spline  interpo¬ 
lation  as  follows: 
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u(x\t)  =  2  0n,tPO){^-  (!) 

where  /3^(j k)= (P\x)^\y)l5^\z)  is  a  separable  cubic 
B  Spline  convolution  kernel,  6n  t  are  the  values  for  the  control 
points  at  phase  t ,  and  Ax  controls  the  grid  spacing  in  each 
dimension. 

To  register  the  CBCT  projections  with  the  planning  CT 
image,  we  adopt  the  metric  of  sum  of  squared  difference, 

$(0)  =  2  -  2  Aijg^2 = 2  -  2  + «)  j2> 

(2) 

where  Yt  is  the  total  attenuation  along  the  ith  projection  ray, 
j  is  the  index  of  the  image  voxel,  and  A  is  the  cone-beam 
projection  matrix.  For  simplicity,  the  phase  index  t  is  ne¬ 
glected  here.  Given  the  projection  data  Y  and  planning  CT 
image  /,  the  deformation  field  u  at  a  particular  phase  can  be 
found  by  minimizing  the  above  objective  function.  General- 
purpose  algorithms  can  be  applied  to  solve  the  optimization 
problem,  and  in  this  work  the  LBFSG  algorithm  is  used,26 
which  requires  the  calculation  of  the  first  derivative  of  the 
cost  function  <F.  It  can  be  shown  that  Eq.  (2)  has  a  closed- 
form  derivative  that  can  be  calculated  explicitly.  For  a  simi¬ 
lar  derivation,  readers  are  referred  to  the  original  work  by 
Zeng  et  al. ,  ’  of  which  the  method  B3  can  be  understood 
as  a  variation. 

II. C.  Simulation  studies 

The  above  model-based  recontruction  method  is  validated 
with  numerical  simulations  under  2D  fan-beam  geometry 
and  3D  cone-beam  geometry.  The  parameters  for  the  simu¬ 
lations  are  similar  to  the  real  setup  of  a  Varian  Trilogy  ma¬ 
chine  (Varian  Medical  Systems,  Palo  Alto,  CA),  where  the 
focal  length  of  the  fan  beam  or  cone  beam  is  150  cm,  and  the 
radius  of  the  detector  rotation  is  50  cm.  The  fan-beam  simu¬ 
lation  has  1024  samples  spaced  by  0.388  mm,  and  2048  pro¬ 
jection  views  were  simulated  with  projection  angles  evenly 
spanned  over  360°,  while  the  cone-beam  simulation  consists 
of  1024X768  pixels  in  each  of  the  total  1200  projections 
with  a  pixel  size  of  0.388  mm.  The  object  used  in  the  fan- 
beam  simulations  came  from  one  slice  of  the  CT  images  of  a 
patient,  which  contained  512  X  512  pixels  with  the  pixel  size 
of  1.0  mm.  For  the  cone-beam  simulation,  the  80  slices  are 
used  with  slicethickness  of  2.5  mm. 

To  include  the  effect  of  the  respiratory-induced  motion, 
we  imposed  an  artificial  perodically  changing  deformation 
field  on  the  image  during  the  simulations.  Specifically,  a 
maximum  deformation  is  determined  by  registering  the  in¬ 
hale  and  exhale  phases  of  a  patient  4DCT  scan.  The  defor¬ 
mation  of  an  intermediate  phase  is  then  obtained  by  linearly 
scaling  down  the  maximum  deformation  field  with  a  scaling 
factor  s  =  cos2( 77^ /  T) ,  where  t  is  the  acquisition  time,  and  T  is 
the  respiration  period.  For  the  fan-beam  case,  the  motion  is 
restricted  in  the  axial  plane,  while  for  the  cone-beam  geom¬ 
etry,  the  motion  is  fully  3D.  The  period  of  the  simulated 
respiratory  motion  is  4  s.  In  reality,  irregular  patient  respira¬ 


tion  is  possible,  which  may  affect  the  accuracy  of  phase  bin¬ 
ning  in  4DCBCT.  However,  this  was  not  considered  in  our 
simulations  because  the  focus  of  this  paper  is  the  reduction 
of  the  view-aliasing  artifacts. 

II.  D.  Patient  study 

A  lung  cancer  patient  had  a  multiple-gantry-rotation 
(MGR)  4DCBCT  scan  in  our  clinic,  and  his  4D  images  were 
used  to  evaluate  our  proposed  approach.  Due  to  the  limit 
access  of  the  available  data,  only  the  method  described  in 
Sec.  2.1.1  was  tested.  The  4DCBCT  was  acquired  with  a 
Trilogy  system  (Varian  Medical  Systems,  Palo  Alto,  CA). 
For  the  MGR  4D  scan,  the  x-ray  tube  current  was  set  to 
32  mA,  and  four  gantry  rotations  were  performed  at  a  speed 
of  6°/s.  Each  rotation  (slightly  over  360°)  consisted  of  over 
680  projections,  and  the  effective  area  of  each  acquired  pro¬ 
jection  image  was  397.312  mm  X  297.984  mm,  containing 
1024  X  768  pixels.  The  average  breathing  cycle  of  the  patient 
was  4.2  s.  The  4DCBCT  images  were  obtained  by  retrospec¬ 
tive  sorting  of  the  MGR  data  into  6  phases  and  subsequently 
reconstructing  the  rebinned  data.  More  details  about  this 
4DCBCT  technique  can  be  found  in  Ref.  16.  To  evaluate  the 
proposed  4D  image  enhancing  method,  the  4D  phases  were 
registered  to  the  0%  phase  with  a  B  Spline  deformable  model 
and  then  superimposed.  The  resultant  phase-0%  image  was 
then  compared  with  the  original  4DCBCT  for  several  axial 
slices  to  demonstrate  the  differences. 

III.  RESULTS 

III.A.  Fan-beam  simulation  study 

The  fan-beam  simulated  projection  data  are  shown  in 
Figs.  2(a)  and  2(b),  where  the  left  is  from  the  static  CT  image 
and  the  right  is  from  the  same  object  with  the  artificial  in¬ 
plane  deformation.  It  is  seen  that  the  object  motion  during 
the  scan  generated  a  large  amount  of  inconsistency  in  the 
projection  data  [Fig.  2(b)].  To  show  its  effect  on  the  resultant 
image,  in  Figs.  2(c)  and  2(d),  the  corresponding  recon¬ 
structed  images  are  compared,  where  artifacts  such  as  blur¬ 
ring,  doubling,  distortion,  and  streak  artifacts  are  observed  in 
the  motion  contaminated  image  [Fig.  2(d)]. 

The  simulated  projections  of  the  continuously  deformed 
image  were  retrospectively  sorted  into  eight  phases.  Five  of 
the  reconstructed  phases  from  the  peak  inspiration  to  the 
peak  expiration  are  illustrated  in  Fig.  3(a),  and  a  zoomed-in 
image  of  phase  1  is  shown  in  Fig.  3(b).  Compared  with  the 
reconstructed  images  in  Figs.  2(c)  and  2(d),  it  can  be  found 
that  the  motion  blurring  is  greatly  reduced,  and  the  boundary 
of  the  primary  tumor  became  much  clearer  after  the  phase 
sorting.  However,  some  low  contrast  structures  (for  example 
the  aorta,  muscle,  etc.)  are  lost  due  to  the  view-aliasing 
streak  artifacts. 

Three  enhanced  images  using  the  proposed  approach  are 
shown  in  Fig.  4  corresponding  to  the  three  motion  model 
derivation  techniques.  The  top  row  in  Fig.  4  shows  the  de¬ 
rived  deformation  between  phase  1  and  phase  4,  where  from 
the  left  to  the  right  are  using  “4D  phase-to-phase  registra- 
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Fig.  2.  Simulated  projections  and  their  corresponding  reconstructed  images.  The  left  column  is  from  the  static  object,  and  the  right  column  is  from  the  same 
object  with  artificial  deformation. 


tion”  (method  Bl),  “CT  image  to  4D  phases”  (method  B2), 
and  “CT  image  to  CBCT  projection”  (method  B3),  respec¬ 
tively.  The  corresponding  reconstructed  images  (at  phase  1) 
are  shown  in  the  bottom  row.  We  find  no  significant  artifacts 
in  the  reconstructed  image,  even  with  certain  errors  in  the 
derived  motion  model  [see  Fig.  4(a)  the  circular  pattern  in 
the  deformation  field  due  to  registering  CBCT  images  con¬ 
taminated  with  view-aliasing  artifacts].  The  image  enhancing 


approach  using  any  of  the  three  motion  models  improved  the 
image  quality  over  the  regular  phase-resolved  image  (Fig.  3), 
and  the  method  as  described  in  B  3  by  matching  CT  image 
with  CBCT  raw  data  generated  the  best  final  images  (less 
streak  artifacts  near  the  boundary  of  the  field  of  view),  be¬ 
cause  a  relatively  more  accurate  motion  model  is  applied. 

In  Figs.  5(a)  and  5(b),  we  quantitatively  compared  the 
CBCT  images  by  plotting  vertical  profiles  passing  through 


Cb) 

Fig.  3.  (a)  Five  examples  of  the  reconstructed  phases  after  the  phase  sorting  of  the  simulated  projection  data  of  the  moving  object;  (b)  a  zoomed-in  image  of 
the  first  reconstructed  phase  where  it  is  found  that  the  motion  blurring  is  reduced  by  phase  binning,  but  some  low-contrast  structures  are  lost  due  to 
view-aliasing  artifacts. 


Medical  Physics,  Vol.  34,  No.  9,  September  2007 


3692 


Li,  Koong,  and  Xing:  Enhanced  4D  CBCT 


3692 


Fig.  4.  Top  row,  from  the  left  to  the  right,  illustrates  the  derived  deformation  field  obtained  by  using  methods  Bl,  B2,  and  B3,  respectively,  and  the  bottom 
row  shows  the  corresponding  reconstructed  images  for  phase  1  using  the  proposed  image  enhancement  technique. 


the  tumor.  It  is  seen  from  Fig.  5(a)  that  the  enhanced 
4DCBCT  (ED_4DCBCT)  images  resulted  from  three  mo- 


Relstive  Position  (mm) 
(a) 


Relative  Position  (mm) 

(b) 

Fig.  5.  Comparison  of  vertical  profiles  passing  through  the  tumor. 


tion  models  have  similar  accuracy  in  terms  of  the  CT  num¬ 
bers.  The  difference  from  the  phantom  (ground  truth)  is 
mainly  due  to  the  filtered  backprojection  reconstruction  pro¬ 
cess.  In  Fig.  5(b),  the  enhanced  4DCBCT  (using  motion 
model  Bl)  is  compared  with  3DCBCT  and  regular  4DCBCT. 
It  can  be  seen  that  the  ED_4DCBCT  results  in  less  noise 
(compared  with  regular  4DCBCT)  while  maintaining  the 
correct  tumor  profile. 

Furthermore,  the  contrast-to-noise  ratio  (CNR)  as  an  im¬ 
age  metric  is  compared  for  the  CBCT  images,  defined  as 
CNR  =  |  (S)  —  (Sq)  |  /  <r,  where  ( S )  and  (S0)  denote  the  average 
signal  within  a  region  of  interest  for  the  tumor  and  soft  tis¬ 
sue,  respectively,  and  a  denotes  the  standard  deviation 
within  the  tumor.  The  calculated  CNR  for  regular  4DCBCT 
is  1.186;  for  enhanced  4DCBCT,  the  CNRs  are  5.487,  5.584, 
and  5.589  for  motion  models  Bl,  B2,  and  B3,  respectively. 

III.B.  Cone-beam  simulation  study 

Figure  6  shows  the  4D  cone-beam  simulation  results  at 
respiratory  phase  1 .  From  the  top  to  the  bottom  rows  are  the 
phantom  images,  the  motion  blurred  3DCBCT  images,  the 
regular  4DCBCT  images,  and  the  enhanced  4DCBCT  im¬ 
ages,  respectively.  From  left  column  to  the  right  column  are 
the  axial,  coronal,  and  sagittal  views.  The  motion  model  ap¬ 
plied  to  the  enhanced  4DCBCT  image  reconstruction  was 
derived  using  method  B2  by  registering  each  phase  of  the 
regular  4DCBCT  to  the  peak-inhale  phantom  images.  Again, 
dramatic  improvement  [Fig.  6(d)  over  Fig.  6(c)]  is  observed 
by  using  the  proposed  approach  for  the  cone-beam  study. 
The  CNRs  for  the  regular  4DCBCT  and  the  enhanced 
4DCBCT  are  0.92  and  4.43,  respectively. 

III.C.  Patient  study 

The  enhanced  4DCBCT  images  were  obtained  for  the  pa¬ 
tient  using  the  proposed  approach  and  the  motion  model  de- 
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Fig.  6.  Simulation  study  using  cone-beam  geometry.  From  left  to  right  are  the  images  of  axial,  coronal,  and  sagittal  views,  respectively.  Row  (a)  shows  the 
4D  phantom;  row  (b)  is  the  3DCBCT  showing  the  motion  blurring  artifacts;  row  (c)  is  the  regular  4DCBCT  images  showing  the  view-aliasing  artifacts;  and 
row  (d)  is  the  enhanced  4DCBCT  images. 


rived  with  method  Bl.  In  Fig.  7  the  resultant  image  at  phase 
1  is  shown  along  with  the  image  of  the  same  phase  before 
the  enhancement.  Because  MGR  technique  was  used  for  the 
4D  CBCT  scan  for  the  patient,  the  view-aliasing  artifacts  in 
the  original  4DCBCT  images  are  not  prominent.  The  im¬ 
provement  by  using  the  proposed  approach  is  not  dramatic  as 
the  previous  simulation  studies.  The  CNRs  are  3.82  and  3.97 
for  the  regular  4DCBCT  and  enhanced  4DCBCT,  respec¬ 
tively.  It  is  found  that  the  proposed  technique  results  in  a 
smoother  image  in  the  uniform  region  leading  to  a  higher 
CNR.  However,  it  is  also  noticed  that  the  process  may  reduce 
the  spatial  resolution.  The  trade-off  between  the  CNR  and 
the  spatial  resolution  using  the  proposed  model-based  recon¬ 
struction  technique  requires  further  systematic  study,  which 
is  beyond  the  scope  of  this  work. 

IV.  DISCUSSION  AND  CONCLUSION 

One  goal  of  4DCBCT  is  to  remove  the  respiratory  motion 
induced  artifacts  in  the  reconstructed  images  to  increase  the 
accuracy  of  target  localization  and  dose  calculation.  Differ¬ 
ent  from  conventional  4DCT, 29-31  4DCBCT  acquisition  with 
an  onboard  imager  usually  has  a  much  slower  speed  of  about 
1  min  per  rotation,  and  due  to  the  limited  total  scan  time  in 
practice,  undersampling  is  often  an  issue  resulting  in  se¬ 
verely  degraded  4D  images,  which  sometimes  are  even 
worse  than  the  motion  contaminated  3D  images.  The  ap¬ 
proach  proposed  in  this  article  reconstructs  the  image  of  any 
particular  phase  by  introducing  additional  information  from 


other  phases,  effectively  increasing  the  sampling  rate  and 
improving  the  image  quality.  In  general,  the  approach  relies 
on  the  accuracy  of  the  motion  model  that  relates  different 
projections.  We  have  demonstrated  with  2D  simulations  that 
the  best  way  to  derive  the  motion  model  is  to  register  the 
artifact- free  CT  images  with  the  CBCT  projections,  if  a  plan¬ 
ning  CT  scan  has  been  performed  prior  to  the  CBCT  scan. 
Though  the  registration  accuracy  in  3D  geometry  may  be 
different  from  2D  geometry,  these  basic  principles  and  rela¬ 
tive  performances  observed  in  this  article  are  expected  true. 
As  an  example,  the  real  patient  case  shown  in  this  article 
demonstrated  the  feasibility  of  the  proposed  method  under 
the  worst  condition  (low  view- aliasing  artifacts  in  the  4D 
images  leave  little  space  for  further  improvement,  and  low 
accuracy  in  motion  derivation  by  using  method  B 1  limits  the 
performance  of  the  proposed  approach). 

As  it  is  known,  the  radiation  dose  has  been  one  of  the 
major  concerns  in  CBCT  imaging  because  of  the  repeated 
use  for  a  given  patient.  With  the  proposed  approach,  it  is 
possible  to  lower  the  x-ray  tube  current,  hence  the  radiation 
exposure  for  the  4DCBCT  scan  while  maintaining  a  mean¬ 
ingful  4D  image,  because  introducing  additional  information 
of  other  phases  during  reconstruction  effectively  increases 
the  photon  statistics.  This  is  similar  to  the  case  of  4D  fan- 
beam  CT  studied  earlier  by  our  group.  The  efficacy  of  the 
approach  in  improving  the  trade-off  between  the  controlled 
radiation  exposure  and  the  resultant  CBCT  image  quality  is 
being  investigated. 
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Fig.  7.  4DCBCT  images  of  a  patient.  The  left  column  shows  three  arbitrary  slices  of  the  peak  inspiration  phase,  and  the  right  column  is  the  corresponding 
images  after  the  proposed  enhancement  with  motion  model  derived  by  using  method  B 1 . 


It  should  be  noted  that  the  proposed  4DCBCT  enhancing 
technique  needs  a  computation  of  the  deformation  field  be¬ 
fore  the  image  reconstruction,  which  takes  extra  time  com¬ 
pared  with  regular  CBCT  imaging.  Currently,  a  3D  image- 
to-image  registration  with  an  image  size  of  512X512X64 
will  take  approximately  30  min.  The  image-to-projection 
registration  method  generally  takes  longer  because  of  the  cal¬ 
culation  of  the  forward  projections.  The  computational  time 
makes  it  difficult  for  online  applications  in  the  treatment 
room.  However,  the  development  of  techniques  in  speeding 
up  the  algorithms  as  well  as  the  computer  hardware  is  highly 
promising,  which  we  believe  will  make  the  4DCBCT  imag¬ 
ing  fast  and  reliable  in  the  near  future.  In  addition,  the 
4DCBCT  image  quality  may  be  further  improved  by  the  ad¬ 
vanced  reconstruction  algorithm  and  helical  scan  mode.33-35 

In  conclusion,  we  have  developed  a  novel  approach  to 
enhance  the  4DCBCT  image  quality.  The  enhancement  is 
achieved  by  increasing  the  angular  sampling  rate  of  each 
phase  of  the  4DCBCT  data  using  information  from  other 
phases,  thus  the  view-aliasing  artifacts  commonly  seen  in  the 
limited-time  scan  4DCBCT  images  are  eliminated  or  signifi¬ 
cantly  reduced.  The  approach  opens  new  opportunities  for 


4DCBCT  to  be  a  more  useful  tool  that  may  substantially 
improve  the  current  cancer  management  in  radiation  oncol¬ 
ogy. 
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OPTIMIZING  4D  CONE-BEAM  CT  ACQUISITION  PROTOCOL  FOR 
EXTERNAL  BEAM  RADIOTHERAPY 
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Purpose:  Four-dimensional  cone-beam  computed  tomography  (4D-CBCT)  imaging  is  sensitive  to  parameters 
such  as  gantry  rotation  speed,  number  of  gantry  rotations,  X-ray  pulse  rate,  and  tube  current,  as  well  as  a 
patient’s  breathing  pattern.  The  aim  of  this  study  is  to  optimize  the  image  acquisition  on  a  patient-specific  basis 
while  minimizing  the  scan  time  and  the  radiation  dose. 

Methods  and  Materials:  More  than  60  sets  of  4D-CBCT  images,  each  with  a  temporal  resolution  of  10  phases, 
were  acquired  using  multiple-gantry  rotation  and  slow-gantry  rotation  techniques.  The  image  quality  was 
quantified  with  a  relative  root  mean-square  error  (RE)  and  correlated  with  various  acquisition  settings; 
specifically,  varying  gantry  rotation  speed,  varying  both  the  rotation  speed  and  the  number  of  rotations,  and 
varying  both  the  rotation  speed  and  tube  current  to  keep  the  radiation  exposure  constant.  These  experiments 
were  repeated  for  three  different  respiratory  periods. 

Results:  With  similar  radiation  dose,  4D-CBCT  images  acquired  with  low  current  and  low  rotation  speed  have 
better  quality  over  images  obtained  with  high  current  and  high  rotation  speed.  In  general,  a  one-rotation 
low-speed  scan  is  superior  to  a  two-rotation  double-speed  scan,  even  though  they  provide  the  same  number  of 
projections.  Furthermore,  it  is  found  that  the  image  quality  behaves  monotonically  with  the  relative  speed  as 
defined  by  the  gantry  rotation  speed  and  the  patient  respiratory  period. 

Conclusions:  The  RE  curves  established  in  this  work  can  be  used  to  predict  the  4D-CBCT  image  quality  before 
a  scan.  This  allows  the  acquisition  protocol  to  be  optimized  individually  to  balance  the  desired  quality  with  the 
associated  scanning  time  and  patient  radiation  dose.  ©  2007  Elsevier  Inc. 

Image-guided  radiotherapy,  Four-dimensional  cone-beam  CT,  Optimization. 


INTRODUCTION 

Onboard  cone-beam  computed  tomography  (CBCT)  imag¬ 
ing  provides  a  convenient  means  for  accurate  patient  setup 
and  dose  verification  (1-8).  However,  when  used  for  thorax 
or  upper  abdomen  imaging,  motion  artifacts  appear  in  the 
reconstructed  images  because  of  intrascanning  organ  mo¬ 
tion  within  the  field  of  view  (FOV).  According  to  the 
International  Electric  Commission  recommendation,  on¬ 
board  CBCT  systems  have  a  limited  gantry  rotation  speed  of 
60  s  per  round.  A  complete  scan  therefore  consists  of 
projection  data  from  10  to  20  respiratory  cycles  of  the 
patient,  resulting  in  large  amount  of  inconsistency  in  the 
CBCT  projection  data.  Reconstruction  algorithms  based  on 
theories  of  static  object  lead  to  images  that  are  significantly 
deteriorated  by  motion  artifacts.  Because  of  the  prolonged 
scanning  time,  these  motion-induced  adverse  effects  in 


CBCT  are  much  more  severe  than  those  seen  in  conven¬ 
tional  CT  scans  (9).  The  artifacts  not  only  blur  the  image, 
but  also  inhibit  direct  use  of  CBCT  for  dose  calculation 
(10-12). 

Four-dimensional  (4D)  CBCT,  or  respiration-correlated 
CBCT,  groups  the  acquired  projection  data  into  several  bins 
as  according  to  their  respiratory  phases.  Each  phase  bin  is 
then  reconstructed  independently  to  obtain  a  volumetric 
image  corresponding  to  that  specific  phase  (13-15).  4D- 
CBCT  not  only  greatly  reduces  the  motion  artifacts  pre¬ 
sented  in  each  of  these  phase-resolved  images,  but  also 
provides  the  dynamic  information  of  the  patient  anatomy, 
which  is  absent  in  the  three-dimensional  (3D)  case,  and 
could  be  very  useful  in  future  4D  radiotherapy  (1). 

Because  the  total  amount  of  projections  acquired  during  a 
4D-CBCT  scan  is  divided  into  several  phase  groups,  the 
number  of  projections  available  for  each  image  reconstruc- 
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tion  becomes  much  less  than  in  the  regular  3D-CBCT  case. 
As  known  from  theory,  the  number  of  projections  needed  to 
avoid  undersampling  artifacts  is  inversely  proportional  to 
the  reconstructed  voxel  size  (16).  An  insufficient  number  of 
projections  will  lead  to  severe  view-aliasing  artifacts.  To 
increase  the  number  of  projections  for  each  phase,  two 
strategies  that  can  be  employed  are:  slowing  down  the 
gantry  rotation  speed  (SGR)  and  multiple  gantry  rotations 
(MGR).  In  addition  to  the  gantry  rotation  speed  and  the 
number  of  gantry  rotations,  the  quality  of  4D-CBCT  images 
is  also  sensitive  to  parameters  such  as  the  X-ray  pulse  rate, 
the  tube  current,  and  the  patient’s  breathing  pattern.  These 
parameters  constitute  a  multidimensional  space,  making 
the  optimal  4D-CBCT  protocol  intractable.  To  balance  the 
tradeoff  between  image  quality,  scan  time,  and  patient  ra¬ 
diation  dose,  a  patient- specific  4D-CBCT  acquisition  pro¬ 
tocol  is  highly  desirable.  This  issue  is  systematically  studied 
in  this  work,  with  the  goal  of  fully  understanding  the  influ¬ 
ences  of  these  various  scanning  parameters  to  optimize 
4D-CBCT  data  acquisition  protocol. 

METHODS  AND  MATERIALS 

4D-CBCT  data  acquisition  system 

An  Acuity  simulator  (Varian  Medical  Systems,  Palo  Alto,  CA) 
was  used  in  this  work  for  all  CBCT  imaging.  In  the  case  of  regular 
3D-CBCT  simulation,  a  gantry  rotation  speed  of  8°/s  is  used.  The 
X-ray  tube  operates  at  125  kVp  and  80  mA  with  a  pulse  width  at 
each  projection  angle  of  9  ms  when  no  bow-tie  filter  is  used.  The 
pulse  width  is  increased  to  25  ms  if  a  bow-tie  filter  is  used  to 
maintain  similar  photon  statistics.  The  data  are  acquired  at  —15 
frames/s  and  a  full  rotation  (slightly  more  than  360°)  consists  of 
about  685  projections  corresponding  to  an  angle  interval  of  about 
0.5°.  The  radiation  dose  of  a  one-rotation  CBCT  scan  at  isocenter 
is  approximately  4.0  cGy.  The  dimension  of  each  acquired  pro¬ 
jection  image  is  397.3  mm  X  298.0  mm,  containing  1,024  X  768 
pixels.  With  a  source-to-axis  distance  of  100  cm  and  source-to- 
detector  distance  of  150  cm,  the  field  of  view  is  around  25  cm  in 
diameter  in  the  transverse  plane  for  a  full-fan  mode.  This  can  be 
close  to  50  cm  if  a  half-fan  mode  is  employed  (by  shifting  the 
flat-panel  detector  laterally).  Here  it  should  be  noted  that  the 
half-fan  mode  generally  results  in  inferior  image  quality  when 
compared  with  that  obtained  with  full-fan  mode  if  all  other  acqui¬ 
sition  parameters  are  maintained  the  same,  and  is  a  consequence  of 
the  loss  of  the  data  redundancy  as  is  in  the  full-fan  mode.  How¬ 
ever,  for  regions  such  as  the  thorax,  a  full-fan  mode  may  result  in 
truncated  images,  which  are  not  appropriate  for  dose  calculation. 
Thus  the  half-fan  mode  is  preferred  when  a  large  field  of  view  is 
required,  and  this  mode  is  therefore  used  in  this  work. 

To  generate  the  4D  image  sets,  the  acquired  CBCT  projections 
must  be  sorted  into  a  number  of  bins  according  to  their  specific 
respiratory  phases.  The  phase  information  of  each  projection  can 
be  obtained  with  the  aid  of  a  motion  tracking  system,  for  example, 
the  Real-time  Position  Management  system  (Varian  Medical  Sys¬ 
tems)  as  done  in  4D-CT  acquisition  with  conventional  diagnostic 
CT  scanners  (17-25),  or  alternatively,  it  can  be  achieved  by 
analyzing  the  acquired  raw  data  in  the  projection  space  (26).  In  the 
latter,  one  or  more  CT-opaque  fiducials  are  placed  on  the  patient 
skin  and  the  locations  of  the  fiducials  in  the  projections  are 
detected  by  a  computer  searching  algorithm.  The  coordinate  of  the 
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Fig.  1.  An  example  of  the  projection- space  phase  tracking  used  in 
this  work.  The  black  solid  squares  and  the  red  open  circles  repre¬ 
sent  the  positions  of  the  fiducial  marker  in  the  real  world  coordi¬ 
nate  system  and  projection  space,  respectively.  The  sinusoidal 
curve  obtained  from  cone-beam  computed  tomography  projections 
can  be  used  to  determine  the  phase  of  each  projection. 


fiducial  is  recorded  as  a  function  of  the  gantry  angle,  revealing  the 
motion  status  of  the  object  for  each  projection.  In  this  work,  the 
projection- space  approach  is  used.  Figure  1  shows  an  example  of 
such  a  plot  for  the  periodically  moving  phantom.  The  projections 
are  phase-tagged  according  to  the  resultant  sinusoidal  curve,  and 
sorted  into  10  phase  groups.  The  data  in  each  phase  bin  are  then 
reconstructed  independently  using  a  Feldkamp  algorithm  with  a 
pixel  size  of  0.5  mm  X  0.5  mm  in  the  cross-section  plane  and  a 
slice  thickness  of  1.0  mm. 

Motion  phantom 

To  investigate  the  influence  of  various  scanning  parameters  on 
the  quality  of  4D-CBCT,  a  motion  phantom  was  constructed.  The 
motion  phantom  consisted  of  a  commercial  CT  calibration  phan¬ 
tom  CatPhan  600  (The  Phantom  Laboratory,  Inc.,  Salem,  NY)  that 
is  placed  on  top  of  a  platform  capable  of  sinusoidal  motion  along 
three  directions:  the  phantom  moved  with  the  maximal  displace¬ 
ments  5.5  cm  in  superoinferior  direction,  1.5  cm  in  anteroposterior 
direction,  and  0.2  cm  in  lateral  direction.  The  period  of  the  motion 
was  continuously  adjustable  in  the  range  of  0.5  s-1  min.  In  this 
work,  three  different  periods  were  used  to  study  the  effects  of 
breathing  cycle  on  the  image  quality,  which  were  3.37  s,  4.59  s, 
and  6.53  s.  4D-CBCT  images  of  the  motion  phantom  were  ac¬ 
quired  using  SGR  or  MGR  methods  as  described  in  the  following 
section.  All  4D  data  acquisitions  were  using  X-ray  tube  current  of 
10  mA  unless  otherwise  stated.  For  comparison,  three  additional 
scans,  with  the  phantom  “frozen”  at  the  “peak-inspiration,”  “mid¬ 
expiration”  and  “peak-expiration”  phases,  were  obtained  as  well 
using  the  standard  3D-CBCT  protocol  (gantry  rotation  speed  at 
8°/s,  X-ray  tube  current  and  voltage  at  80  mA  and  125  kVp, 
respectively).  The  three  scans  served  as  our  “standards”  for  these 
phases. 

4D-CBCT  image  analysis 

Because  the  images  in  this  work  involve  significant  artifacts,  a 
common  metric  such  as  contrast- to-noise  ratio  (CNR),  defined  as 
CNR  =  I  <  S{  —  S2  >  l/cr,  may  not  be  suitable,  in  the  sense  that 
a  high  CNR  output  may  not  represent  a  better  image.  For  example, 
dark  streak  artifacts  may  accidentally  decrease  the  mean  value  of 
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S2,  leading  to  a  high  CNR  value.  Thus  the  CNR  metric  may  highly 
depend  on  the  selected  regions  of  interest.  To  have  a  more  robust 
and  accurate  assessment  of  the  images  obtained  using  different 
CBCT  settings,  a  “relative  root  mean  square  error”  (denoted  by 
RE)  is  chosen  as  the  figure  of  merit  of  the  image  quality,  which  is 
defined  as 


where  $4D  denotes  the  4D  single-phase  image,  and  S0  is  the 
standard  80-mA  3D-CBCT  images  of  the  phantom  frozen  at  the 
same  phase.  The  summation  runs  over  all  voxels  of  the  images.  In 
Eq.  1,  S0  image  is  used  as  the  “gold  standard,”  and  the  mean  square 
error  between  the  4D  images  and  the  gold  standard  is  normalized 
to  the  mean  square  of  the  true  intensity. 

Slow-gantry  rotation  strategy 

A  sufficient  number  of  projections  must  be  collected  for  each 
phase  to  remove  the  view-aliasing  artifacts  in  4D-CBCT.  One  way 
to  achieve  this  is  to  slow  down  the  gantry  rotation.  In  this  way,  the 
number  of  breathing  cycles  covered  by  a  full  rotation  is  increased, 
leading  to  more  projections  in  each  phase  group.  Figure  2  is  a 
sketch  illustrating  the  relationship  between  the  gantry  rotation 
schemes  and  the  breathing.  The  dark  areas  represent  the  available 
projections  that  belong  to  the  same  phase  group. 

To  examine  the  influence  of  gantry  rotation  speed  on  the  result¬ 
ant  image  quality,  scans  were  made  of  the  motion  phantom  for 
eight  gantry  rotation  speeds  ranging  from  87s  to  17s.  All  other 
system  variables  were  kept  constant  with  the  X-ray  tube  current 
and  voltage  in  these  eight  scans  being  held  at  10  mA  and  125  kVp, 
respectively.  This  increased  the  total  number  of  projections  from 
about  680  for  87s  speed  to  about  5,440  for  17s  speed.  Conse¬ 
quently,  the  number  of  projections  in  each  phase  bin  increased  by 
~8  times.  The  resultant  4D-CBCT  images  were  then  analyzed 
using  the  RE  metric  as  defined  in  Eq.  1. 

Another  important  quantity  for  SGR  4D  acquisition  is  the  rela¬ 
tion  between  image  quality  and  the  associated  radiation  exposure. 
For  a  constant  X-ray  pulse  width,  the  radiation  dose  is  directly 
proportional  to  the  number  of  projections  (beam-on  time)  and 
essentially  linear  to  the  X-ray  tube  current.  Hence  a  scan  with  a 
rotational  speed  of  17s  at  10  mA  current  will  result  in  a  similar 
amount  of  radiation  as  a  scan  with  a  27s  rotation  speed  at  doubled 
the  tube  current.  To  investigate  this  image  quality-tube  current 
relationship,  4D-CBCT  images  of  the  motion  phantom  were  ac¬ 
quired  for  17s- 10  mA,  27s-20  mA,  47s-40  mA,  57s-50  mA,  and 
87s-80  mA,  for  each  of  the  three  motion  periods  mentioned 
previously.  These  scanning  schemes  deliver  a  similar  radiation 
dose  to  the  phantom.  A  quantitative  comparison  of  the  4D  images 
was  carried  out. 

Multiple- gantry  rotation  strategy 

An  alternative  to  SGR  is  MGR  4D  acquisition,  in  which  the 
gantry  now  rotates  multiple  times  back  and  forth.  As  with  SGR, 
the  projections  of  the  same  phase  are  again  collected  and  used  for 
4D-CBCT  image  reconstruction.  The  MGR  acquisition  scheme  in 
relation  to  the  breathing  pattern  is  sketched  in  Fig.  2c.  As  shown, 
both  SGR  and  MGR  strategies  lead  to  an  increase  in  the  number  of 
projections  over  conventional  CBCT  (Figs.  2b,  2c);  however,  the 
angular  distribution  of  the  projections  of  a  given  phase  group  is 
very  different.  In  SGR,  the  projections  belonging  to  a  particular 


breathing  phase  are  regularly  distributed  in  the  angular  space, 
whereas  in  MGR,  two  or  more  rotations  lead  to  the  projection 
becoming  irregular.  As  a  result,  now  there  is  a  chance  that,  during 
the  multiple  gantry  rotations,  projections  at  the  same  phase  and 
gantry  angle  may  occur  more  than  once,  leading  to  redundancy  of 
the  data  and  resulting  in  less  effective  4D-CBCT  reconstruction. 
The  angular  overlap  or  partial  overlap  of  the  projection  data  can, 
in  principle,  be  avoided/reduced  by  properly  selecting  the  initial 
phase  of  the  two  (or  multiple)  scans.  Generally,  as  illustrated  in 
Fig.  2d,  the  chance  for  two  successive  scans  to  overlap  is  minimal 
when  the  gantry  moves  in  opposite  directions.  One  can  imagine 
that  if  two  scans  are  in  the  same  direction,  at  every  angle,  the  two 
projections  from  the  two  scans  may  be  at  same  phase  in  the 
extreme  case.  This  scenario  will  never  happen  if  the  two  scans  are 
in  opposite  directions.  Two-rotation  MGR  scans  with  the  gantry 
“back  and  forth”  were  performed  for  each  of  the  previously  men¬ 
tioned  three  motion  periods  and  the  results  were  compared  with 
that  obtained  using  SGR.  To  maintain  a  constant  radiation  expo¬ 
sure,  the  gantry  rotation  speed  in  the  MGR  approach  was  increased 
by  a  factor  of  two  as  compared  with  that  of  a  single  rotation  scan. 

Relative  gantry  speed 

Four-dimensional-CBCT  image  quality  depends  not  only  on  the 
gantry  rotation  speed,  but  also  on  the  patient  respiratory  pattern. 
This  suggests  that  a  patient-specific  setting  might  be  needed  to 
achieve  a  similar  quality  of  4D-CBCT  images.  Here  we  propose  a 
new  measure  of  relative  speed  to  combine  the  two  variables  of 
gantry  rotation  and  patient  respiration,  which  is  defined  as  follows 

Relative  Speed  =  Respiratory  Period(s) 

X  Gantry  Rotation  Speed(deg/s).  (2) 

Relative  speed  has  a  unit  of  degree,  which  is  the  spanned  angle 
of  the  CBCT  scan  during  the  time  of  one  respiratory  cycle. 

The  data  acquired  in  these  studies,  for  different  motion  periods 
and  different  scanning  speed  at  constant  X-ray  tube  current  (10 
mA)  and  voltage  (125  kVp),  were  reanalyzed  to  illustrate  the 
utility  of  the  new  concept  in  optimizing  the  4D-CBCT  image 
acquisition  protocol  on  a  patient- specific  basis. 

RESULTS 

4D-CBCT  image  quality  as  a  function  of  gantry 
rotation  speed 

Four-dimensional-CBCT  images  acquired  with  eight  dif¬ 
ferent  gantry  rotation  speeds  and  constant  X-ray  tube  cur¬ 
rent  (10  mA)  for  the  phantom  moving  at  a  period  of  3.37  s 
are  summarized  in  Fig.  3.  From  the  top-left  to  the  bottom- 
middle,  the  images  correspond  to  gantry  rotation  speeds 
from  87s  to  17s.  Without  loss  of  generality,  only  the  phase 
0%  images  of  the  4D-CBCT  acquisitions  are  compared 
here.  Note  that  all  other  imaging  parameters  such  as  the  tube 
current  and  voltage  were  the  same  for  all  the  eight  scans. 
The  golden  standard  image  (obtained  with  a  standard  3D- 
CBCT  protocol  and  without  motion)  is  also  shown  at  the 
bottom  right  in  Fig.  3.  As  expected,  it  is  found  that  the 
image  quality  improves  with  decreasing  gantry  rotation 
speeds.  Similar  trends  were  also  observed  for  the  phantom 
moving  with  the  other  two  periods  (4.93  s  and  6.53  s). 
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Fig.  2.  Illustration  of  four-dimensional  cone-beam  computed  tomography  acquisitions  with  gantry  rotation  speed  (SGR) 
and  multiple  gantry  rotations  (MGR)  methods.  The  shadow  in  the  diagrams  represents  the  projections  available  for 
reconstruction  of  one  particular  phase,  (a)  Single  rotation  with  rotation  speed  of  a  conventional  three-dimensional  scan, 
in  which  the  number  of  projections  in  a  phase  bin  is  relatively  small;  (b)  with  SGR,  a  scan  expands  more  respiratory 
cycles,  hence  more  projections  of  the  same  phase  can  be  collected;  (c)  with  MGR,  the  number  of  projections  in  a  phase 
group  increases  as  shown  by  the  gray  and  blue  shadow;  (d)  the  phase  distribution  in  two  successive  rotations  of  the  MGR 
method.  Note  that  if  two  projections  overlap  at  a  particular  angle,  for  example,  for  0%  phase,  the  projections  at  neighbor 
angles  will  be  in  different  phases  because  the  two  rotations  are  in  opposite  directions. 


The  quantitative  relation  between  image  quality  and  gan¬ 
try  rotation  speed  was  analyzed  using  the  image  metric  of 
REs  and  the  results  were  plotted  in  Fig.  4.  For  all  three 


respiratory  periods,  the  relative  errors  are  found  to  increase 
with  increasing  gantry  rotation  speeds,  indicating  a  drop  in 
the  image  quality.  Also  from  Fig.  4,  it  is  seen  that  at  a  given 
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Fig.  3.  From  the  top-left  to  the  bottom-right  are  the  four-dimen¬ 
sional  cone-beam  computed  tomography  (4D-CBCT)  images  of 
phase  0%  obtained  by  varying  the  gantry  rotation  speed  from  87s 
to  17s,  and  the  three-dimensional  CBCT  images  of  the  phantom 
“frozen”  at  the  same  phase.  All  4D-CBCT  were  acquired  with 
X-ray  tube  current  of  10  mA. 


gantry  rotation  speed,  the  image  quality  improves  (smaller 
REs)  with  increasing  phantom  velocities.  This  is  not  sur¬ 
prising  because  a  full  scan  at  fast  phantom  speed  will 
contain  more  respiratory  cycles,  leading  to  a  reduction  of 
the  gap  (Fig.  2)  between  the  projections  of  same  phase  in 
two  successive  cycles,  hence  less  view-aliasing  artifacts.  To 
help  visualize  this  effect,  Fig.  5  shows  the  0%  phase  images 
obtained  for  the  three  different  motion  modes  with  the  same 
gantry  rotation  speed  of  27s. 
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Fig.  4.  Relation  between  relative  root  mean- square  error  and 
gantry  rotation  speed  for  three  motion  modes  with  periods  of 
3.37  s,  4.93  s,  and  6.53  s,  respectively. 


Fig.  5.  Four-dimensional  cone-beam  computed  tomography  im¬ 
ages  of  phase  0%  obtained  at  the  same  gantry  rotation  speed  of 
27s,  but  different  phantom  motion  periods. 


SGR  acquisition  under  the  condition  of  a  constant 
radiation  exposure 

As  mentioned  previously,  the  photon  flux,  noise  level  and 
patient  radiation  exposure  are  all  related  to  the  tube  current. 
This  is  investigated  by  comparison  of  the  0%  phase  4D 
images  for  scans  of  87s  at  80  mA,  57s  at  50  mA,  47s  at  40 
mA,  27s  at  20  mA,  and  17s  at  10  mA.  Because  the 
radiation  exposure  is  approximately  proportional  to  the  tube 
current,  scans  with  these  eight  settings  deliver  a  similar 
level  of  radiation  dose  to  the  phantom.  Again,  the  image 
metric  of  RE  was  employed  to  quantify  the  difference 
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Fig.  6.  The  effects  of  increasing  tube  current  and  gantry  rotation 
speed.  With  the  same  radiation  dose,  lower  speed  results  in  higher 
quality  four-dimensional  images  as  a  result  of  increased  number  of 
projections  per  phase. 


among  the  images.  As  shown  in  Fig.  6,  the  image  acquired 
with  the  combination  of  l°/s  and  10  mA  setting  had  the  best 
quality  (smallest  RE).  For  visual  comparison,  the  l°/s — 10 
mA  and  8°/s — 80  mA  images  are  shown  in  Fig.  7.  Moving 
to  slower  gantry  rotation  speed  will  lead  to  each  phase 
accumulating  more  projection  data  that  can  be  used  for 
reconstruction  as  compared  with  faster  scans.  As  a  conse¬ 
quence,  the  integral  photon  number  (filtered  back  projec¬ 
tion)  at  a  spatial  point  is  not  reduced  in  the  SGR  low-current 
acquisition  as  compared  with  high-current,  high-speed 
scans.  Furthermore,  at  low  scan  speed,  the  data  become 
more  evenly  spared  along  the  360°  circle  for  each  phase, 
leading  to  less  view-aliasing  artifact  and  consequently  to 
higher  image  quality. 

As  shown  in  Fig.  6,  it  is  found  that  RE  increases  almost 
linearly  as  a  function  of  the  scan  speed.  Because  the  total 
number  of  projections  in  a  scan,  and  therefore  the  number  of 
projections  in  each  phase,  varies  linearly  with  the  gantry 
rotation  speed,  it  seems  that  the  RE  measure  is  an  appro¬ 
priate  image  metric,  which  has  a  linear  relationship  to  the 
number  of  projections  available  for  generating  the  phase 
image. 

4D-CBCT  imaging  with  single-  and 
double-gantry  rotations 

Both  SGR  and  MGR  resulted  in  an  increased  number  of 
projections.  However,  as  illustrated  in  Figs.  2b  and  2c,  the 
angular  distribution  of  the  projections  is  different  for  the  two 
acquisition  schemes.  In  general,  the  efficiency  of  the  MGR 
method  depends  on  the  level  of  projection  redundancy  or 
overlap  of  the  data  from  the  two  or  more  gantry  rotations. 
The  4D-CBCT  images  acquired  with  SGR  and  MGR  were 
compared  using  the  RE  metric,  and  the  results  are  summa¬ 
rized  in  Fig.  8.  Here  it  is  found  that  for  all  three  motion 
modes  (3.37  s,  4.93  s,  and  6.53  s,  respectively)  and  for 
different  combinations  of  the  gantry  rotations  speed  and 
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number  of  rotations,  SGR  almost  always  led  to  a  lower  RE 
or  better  image  quality  than  the  MGR  method. 

4D-CBCT  image  quality  as  a  function  of  the 
relative  speed 

The  relationship  between  the  image  quality  and  the  rela¬ 
tive  speed  is  presented  in  Fig.  9.  The  data  shown  in  Fig.  9 
were  fitted  with  a  second  order  polynomial.  The  parameters 
of  the  model  and  the  goodness  of  the  fit  are  also  shown  in 
the  figure.  The  monotonic  dependence  of  the  RE  on  the 
relative  speed  suggests  that  the  relative  speed  is  a  meaning¬ 
ful  concept  in  characterizing  the  4D-CBCT  data  acquisi¬ 
tion.  The  quantity  “condenses”  two  important  parameters 
(patient-specific  respiration  period  and  the  gantry  rotation 
speed)  into  one  and  is  useful  for  us  to  optimize  the  4D- 
CBCT  acquisition  protocol.  In  reality,  because  the  respira¬ 
tory  period  of  a  patient  can  be  measured  with  the  aid  of  a 
real-time  position  management  signal  or  similar  device 
( e.g .,  strain  gauge  signal)  before  the  4D-CBCT  scanning,  a 
suitable  gantry  speed  can  be  derived  to  achieve  a  prespeci¬ 
fied  image  quality  before  4D-CBCT  scan.  Alternatively,  one 
can  use  the  curve  as  shown  in  Fig.  9  to  preestimate  the  4D 
image  quality  for  an  acquisition  setting  and  to  estimate  the 
corresponding  patient  radiation  dose. 

DISCUSSION 

Cone-beam  CT  volumetric  imaging  integrated  with  a 
medical  linear  accelerator  opens  new  avenues  for  improving 
current  radiation  oncology  practice.  In  reality,  CBCT  has 
two  important  applications:  patient  setup  and  dose  recon¬ 
struction/verification.  By  imaging  the  patient  routinely  dur¬ 
ing  a  course  of  radiation  therapy,  the  accuracy  of  the  patient 
setup  can  potentially  be  improved.  Furthermore,  the  CBCT 
provides  a  pretreatment  patient  model  on  which  the  dose 
calculation  can  be  performed  using  intended  fluence  maps 
from  the  planning  system  or  other  means.  Both  applications 
rely  on  the  fidelity  and  quality  of  the  volumetric  images. 
When  imaging  in  the  region  of  thorax  and  upper  abdomen, 
however,  respiration-induced  artifacts,  such  as  blurring, 
doubling,  streaking,  and  distortion  are  observed,  which 
heavily  degrade  the  image  quality,  and  affect  the  target 


Fig.  7.  Comparison  of  four-dimensional  cone-beam  computed  to¬ 
mography  images  acquired  with  speed  l°/s  tube  current  10  mA 
(left)  and  speed  87s  tube  current  80  mA  (right). 


Optimizing  4D  cone-beam  CT  •  T.  Li  and  L.  Xing 


1217 


Fig.  8.  The  effects  of  increasing  the  number  of  rotations  and  gantry  speed.  With  the  same  radiation  dose,  for  most  of 
the  cases,  low  rotation  speed  results  in  a  better  image  quality  than  multiple  rotations. 


localization  ability  and  the  accuracy  of  dose  verification. 
These  artifacts  in  CBCT  are  much  more  severe  than  those 
found  in  conventional  CT  exams  because  of  the  relatively 
slow  gantry  rotation.  In  conventional  CT,  each  rotation  of 
the  scan  can  be  completed  within  0.5  s,  during  this  period 
the  organ/tumor  motion  is  relatively  small.  Patient  body 
restraints  and  breath-hold  techniques  can  be  used  to  mini¬ 
mize  the  motion  if  necessary.  On  the  contrary,  in  CBCT 
scan,  the  gantry  rotation  speed  is  much  slower,  typically 
around  1  min  for  a  full  360°  scan  in  acquiring  the  projection 
data,  which  covers  more  than  10  breathing  cycles  for  most 
patients.  4D-CBCT  or  phase-correlated  CBCT  is  an  effec¬ 
tive  way  to  reduce/eliminate  the  motion  artifacts  and  makes 
it  useful  for  guiding  the  patient  setup  and  dose  validation  in 
the  presence  of  organ  motion. 

Radiation  dose  is  an  important  issue  in  CT  imaging  and, 
in  particular  the  onboard  CBCT  imaging  because  of  the 
repeated  use  of  modality  for  a  given  patient.  We  found  that, 
as  we  increase  the  number  of  projections  per  phase  by 
slowing  down  the  gantry  rotation  speed  or  multiple  gantry 
rotations,  the  tube  current  can  be  lowered  accordingly.  In 
this  way,  the  4D  image  quality  is  generally  not  compro¬ 


mised  and  one  can  obtain  decent  4D-CT  images  without 
increasing  the  radiation  exposure.  In  a  sense,  this  scheme  is 
to  “spread”  the  photons  of  a  projection  in  3D-CBCT  to  a 
range  of  angles  in  4D-CBCT  imaging,  which  represents  a 
better  tradeoff  between  the  signal-to-noise  ratio  and  the 
reduction  in  motion  artifacts.  In  addition  to  the  reduced 
patient  dose,  low  current  is  essential  to  avoid  overheating 
the  X-ray  tube,  especially  for  an  onboard  imager  without  oil 
cooling  system.  Although  the  4D  image  quality  directly 
relates  to  the  number  of  projections  available  for  each 
phase,  there  is  a  limit  for  slowing  down  gantry  speed. 
Beyond  this  limit,  the  time  needed  to  complete  the  scan  may 
be  too  long  to  cause  patient  discomfort  and  extra  intrascan 
motion  (other  than  the  respiration-induced  motion)  may 
lead  to  undesirable  artifacts.  There  is  a  practical  need  to 
keep  the  scanning  time  as  short  as  possible  while  maintain¬ 
ing  an  acceptable  image  quality.  The  result  shown  in  Fig.  9 
provides  guidance  for  the  4D-CBCT  acquisition  on  a  pa¬ 
tient-specific  basis,  which  predicts  the  quality  of  the  4D 
images  for  a  certain  gantry  rotation  speed  and  patient  respi¬ 
ratory  period.  Though  patient  irregular  breathing  patterns 
may  alter  the  number  and  angular  distribution  of  the  pro- 
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Fig.  9.  Relationship  between  image  quality  and  relative  speed  of 
the  four-dimensional  cone-beam  computed  tomography  scan.  The 
data  were  acquired  for  three  motion  modes  with  gantry  rotation 
speed  varying  from  l°/s  to  8°/s  and  constant  X-ray  tube  current 
(10  mA)  and  voltage  (125  kVp). 

jections  in  each  phase  and  thus  change  the  image  qualities, 
Fig.  9  will  nonetheless  represent  an  approximate  relation 
between  the  image  quality  and  the  average  patient  respira¬ 
tory  period.  Obviously,  this  approximation  is  valid  only 
when  the  irregularity  of  the  breathing  pattern  is  not  far  from 
ideal  periodic  motion.  It  is  also  worth  noting  that  the  influ¬ 
ence  of  breathing  irregularity  on  the  4D  images  is  different 
from  that  in  conventional  fan  beam  CT,  in  which  the  bin¬ 
ning  artifacts  may  arise  due  to  mismatch  of  the  slices 
corresponding  to  different  couch  positions. 

Although  4D-CBCT  technique  (SGR  or  MGR)  with 
phase  binning  before  reconstruction  reduces  the  motion 
artifacts  and  improves  the  image  quality,  it  deals  with  each 
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phase  independently  and  ignores  the  correlation  of  the  pro¬ 
jections  among  different  phase  groups.  It  is  of  great  interest 
to  reconstruct  a  4D  object  (with  spatial  and  temporal  pa¬ 
rameters)  rather  than  a  series  of  independent  3D  objects 
(with  only  spatial  parameters),  by  using  simultaneously  all 
the  projection  data  collected  in  different  respiratory  phase 
groups.  With  the  aid  of  deformable  registration  and  novel 
reconstruction  techniques  (9,  22,  27-33),  it  is  possible  to 
incorporate  all  projection  data  of  different  phases  into  one 
image  reconstruction  process.  In  this  way,  information  from 
the  measurement  (CBCT  scan)  is  fully  used  in  the  data 
processing,  and  more  accurate  imaging  and  better  dose 
efficiency  may  be  achieved. 

CONCLUSIONS 

We  have  demonstrated  that  4D-CBCT  images  with  sig¬ 
nificantly  reduced  artifacts  can  be  obtained  with  a  commer¬ 
cially  available  on-board  imager  system.  A  measure  of 
relative  error  was  used  to  quantify  the  image  quality.  At  the 
same  radiation  dose,  the  low-current,  low-speed  acquisition 
is  superior  to  acquisition  with  a  high  current  and  high  speed 
(<8°/s),  and,  in  most  cases,  for  the  same  total  scanning 
time,  one-rotation  low-speed  scan  is  more  preferable  as 
compared  with  an  acquisition  based  on  multiple  rotations 
and  higher  speed.  Furthermore,  we  suggested  that  4D- 
CBCT  gantry  rotation  speed  should  be  determined  on  a 
patient- specific  basis  and  have  established  a  quantitative 
relation  between  the  4D-CBCT  image  quality  and  the  rela¬ 
tive  speed.  The  approach  optimizes  the  4D  acquisition  and 
makes  an  efficient  use  of  the  available  4D-CBCT  technol¬ 
ogy.  The  study  lays  foundation  for  clinical  application  of 
4D-CBCT  and  should  have  significant  implication  to  image- 
guided  radiation  therapy. 
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The  purpose  of  this  work  is  to  develop  a  novel  strategy  to  automatically  map  organ  contours  from 
one  phase  of  respiration  to  all  other  phases  on  a  four-dimensional  computed  tomography  (4D  CT). 

A  region  of  interest  (ROI)  was  manually  delineated  by  a  physician  on  one  phase  specific  image  set 
of  a  4D  CT.  A  number  of  cubic  control  volumes  of  the  size  of  ~1  cm  were  automatically  placed 
along  the  contours.  The  control  volumes  were  then  collectively  mapped  to  the  next  phase  using  a 
rigid  transformation.  To  accommodate  organ  deformation,  a  model-based  adaptation  of  the  control 
volume  positions  was  followed  after  the  rigid  mapping  procedure.  This  further  adjustment  of 
control  volume  positions  was  performed  by  minimizing  an  energy  function  which  balances  the 
tendency  for  the  control  volumes  to  move  to  their  correspondences  with  the  desire  to  maintain 
similar  image  features  and  shape  integrity  of  the  contour.  The  mapped  ROI  surface  was  then 
constructed  based  on  the  central  positions  of  the  control  volumes  using  a  triangulated  surface 
construction  technique.  The  proposed  technique  was  assessed  using  a  digital  phantom  and  4D  CT 
images  of  three  lung  patients.  Our  digital  phantom  study  data  indicated  that  a  spatial  accuracy  better 
than  2.5  mm  is  achievable  using  the  proposed  technique.  The  patient  study  showed  a  similar  level 
of  accuracy.  In  addition,  the  computational  speed  of  our  algorithm  was  significantly  improved  as 
compared  with  a  conventional  deformable  registration-based  contour  mapping  technique.  The  ro¬ 
bustness  and  accuracy  of  this  approach  make  it  a  valuable  tool  for  the  efficient  use  of  the  available 
spatial-tempo  information  for  4D  simulation  and  treatment.  ©  2007  American  Association  of 
Physicists  in  Medicine.  [DOI:  10.1118/1.2780105] 

Key  words:  CT,  image  registration,  contour  mapping,  IGRT 


I.  INTRODUCTION 

A  longstanding  question  in  radiation  therapy  is  how  to  accu¬ 
rately  and  efficiently  segment  a  region  of  interest  (ROI)  such 
as  a  tumor  target  volume  or  a  critical  structure.1'6  In  spite  of 
intense  research  efforts  in  the  past  few  decades,  ROI  seg¬ 
mentation  remains  a  time  consuming  task  in  treatment  plan¬ 
ning.  In  most  cases,  the  segmentation  is  performed  manually 
in  a  slice-by-slice  fashion,  creating  a  strong  need  for  auto¬ 
mated  segmentation  tools  in  the  clinics.  The  introduction  of 
four-dimensional  computed  tomography  (4D  CT)  in  radia¬ 
tion  oncology  practice  further  amplifies  this  need  as  the  num¬ 
ber  of  images  to  be  segmented  is  increased  dramatically.7"15 
Generally,  a  4D  CT  scan  consists  of  5-10  sets  of  three- 
dimensional  (3D)  CT  images,  each  representing  the  patient 
anatomy  at  a  specific  phase  of  respiration.  For  4D  radiation 
therapy  applications,  it  is  labor  intensive  to  follow  the  3D 
approach  of  manual  segmentation  due  to  the  immense  work¬ 
load  associated  with  this  process. 

A  natural  way  to  deal  with  the  4D  segmentation  problem 
is  to  start  with  a  known  set  of  contours  for  a  selected  phase 
and  map  these  contours  onto  all  other  phases.  ROIs  for  the 
selected  phase  are  first  manually  contoured,  similar  to  that 
done  in  treatment  planning  based  on  3D  CT  image  data  sets. 
The  mapping  procedure  can  be  accomplished  with  the  aid  of 
a  computer  algorithm  that  registers  an  arbitrary  point  on  the 
selected  phase  to  the  corresponding  points  on  all  other 


phases.  While  conceptually  simple,  the  implementation  of 
this  idea  is  not  straightforward.  An  intelligent  algorithm  ca¬ 
pable  of  providing  accurate  point-to-point  correspondence 
between  the  phased  images,  or  at  least  between  points  within 
the  ROIs,  is  the  key  to  the  success  of  this  approach.  Various 
studies  have  investigated  algorithms  to  automatically  map 
contours  using  deformable  image  registration  and  surface 
mapping  techniques.  One  method  is  based  on  a  deformable 
registration  model.16"18  This  method  has  limited  accuracy, 
especially  in  the  regions  proximate  to  the  interfaces  of  dif¬ 
ferent  organs,  and  is  brute-force  in  nature,  which  entails  a 
large  amount  of  computations.  In  reality,  contour  mapping  is 
a  regional  problem  and  a  global  association  of  the  phase- 
based  images  is  neither  necessary  nor  efficient.  Surface  map¬ 
ping  achieves  the  stated  goal  of  contour  transformation  by 
iteratively  deforming  the  ROI  contour-extended  surface  until 
the  optimal  match  with  the  reference  is  reached.2,19'22  Nu¬ 
merous  surface  mapping  techniques,  such  as  spatial  parti¬ 
tioning,  principal  component  analysis,  conformal  mapping, 
rigid  affine  transformation,  deformable  contours,  and  warp¬ 
ing  based  on  the  thin-plate  spline  (TPS),  have  been  devel¬ 
oped  over  the  years  and  the  end  point  of  all  of  these  tech¬ 
niques  is  a  mapping  between  topological  components  of  the 
input  surfaces  that  allow  for  transfer  of  annotations.  This 
type  of  computation  is  inherently  more  efficient  in  compari¬ 
son  with  the  deformable  model-based  approaches,  but  it  suf- 
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fers  from  the  fact  that  the  resultant  mapping  heavily  depends 
on  the  model  used  and  the  fact  that  the  model  parameters  in 
the  calculations  are  not  physically  transparent. 

In  this  work,  we  combine  the  useful  features  of  the  two 
different  types  of  techniques  for  contour  mapping.  Our  work¬ 
ing  hypothesis  is  that  information  contained  in  the  boundary 
region  is  often  sufficient  to  guide  the  contour  mapping  pro¬ 
cess  without  relying  on  the  use  of  an  ad  hoc  surface  deform¬ 
ing  model.  In  the  proposed  technique,  the  neighborhood  in¬ 
formation  of  the  contour  surface  is  captured  by  a  series  of 
small  cubic  (or  other  shaped)  control  volumes  placed  around 
the  surface.23  The  collection  of  the  centers  of  the  control 
volumes  is  a  representation  of  the  contour  surface.  The  ROI 
contour  mapping  proceeds  iteratively  under  the  guidance  of 
the  information  contained  in  the  control  volumes.  The  pro¬ 
posed  method  is  illustrated  by  a  digital  phantom  experiment 
and  three  clinical  lung  case  studies.  The  results  suggest  that 
the  technique  is  capable  of  automatically  mapping  contours 
among  the4D  CT  phases  with  clinically  acceptable  accuracy. 


Fig.  1.  A  schematic  drawing  of  the  placement  of  control  volumes  (shown  as 
squares)  on  the  manually  segmented  contour  (depicted  curve).  Each  control 
volume  is  typically  ~1  cm  in  cubic  shape  as  shown  in  the  zoomed  image. 
No  interpolated  volumes  are  displayed. 


II.  MATERIALS  AND  METHOD 
II.A.  Software  platform 

The  Insight  Toolkit24  (ITK)  and  the  Visualization 
Toolkit25  (VTK),  which  are  open  source  cross-platform 
C  +  +  software  toolkits  and  are  freely  available  for  research 
purposes,  were  used  in  this  study.  A  variety  of  methods  have 
been  programmed  into  the  ITK  platform  for  image  registra¬ 
tion  and  segmentation.  ITK  was  used  for  automatic  mapping 
while  VTK  was  mainly  used  for  3D  visualization  and  con¬ 
tour  construction. 

II.B.  Image  acquisition 

The  4D  CT  image  data  sets  for  three  lung  cancer  patients 
were  acquired  with  a  multislice  helical  CT  scanner  (Discov¬ 
ery  ST,  GE  M  edical  System,  M  ilwaukee,  Wl).  The  collected 
data  were  sorted  into  ten  phase  bins.10  The  4D  CT  image  sets 
for  all  patient  studies  were  reconstructed  with  a  2.5  mm  slice 
thickness.  The  size  and  pixel  resolution  of  each  CT  slice  was 
512x512  and  0.98x0.98  mm2,  respectively.  The  4D  CT 
images  were  transferred  through  DICOM  to  a  personal  com¬ 
puter  with  a  Pentium  IV  (2.66  GHz)  processor  for  image 
processing.  One  of  the  phases,  which  is  referred  to  as  the 
template  phase,  was  selected  for  manual  segmentation  of  the 
ROIs.  We  refer  to  the  other  phases  in  the  4D  CT  image  set  as 
the  target  phases  with  corresponding  target  contours.  The 
enrolled  patients  in  this  study  were  under  an  Institutional 
Review  Board  approved  protocol. 

II.C.  Placement  of  control  volumes  along  the  ROI 
contour  surface 

The  task  of  mapping  a  contour  is  to  find  its  corresponding 
location  on  the  target  phase  for  an  arbitrary  point  on  the 
contour  drawn  on  the  template  image.  In  general,  the  image 
feature  surrounding  a  contour  point  can  be  used  as  a  signa¬ 
ture  of  the  point  to  aid  the  search  for  its  corresponding  po¬ 
sition  on  the  target  phase.  In  this  study,  the  image  feature  at 


each  point  was  captured  by  introducing  a  cubic  control  vol¬ 
ume  (~1  cm  in  size)  centered  at  the  point.  The  manual  ROI 
delineation  was  first  performed  on  the  template  phase  using 
theVarian  Eclipse  TPS  (Varian  Medical  Systems,  Palo  Alto, 
CA).  Afterwards,  the  contours  were  exported  from  the  TPS 
to  a  local  computer  for  contour  propagation.  The  exported 
contours  are  polygons  on  individual  slices,  and  the  vertices 
of  the  polygons  were  used  as  the  locations  of  control  vol¬ 
umes.  This  is  depicted  in  Fig.  1,  where  the  curve  represents 
the  contour  and  the  squares  represent  the  control  volumes. 
The  center  of  each  control  volume  was  set  at  a  contour  point, 
typically  0.5-1  cm  from  the  next  point.  More  control  vol¬ 
umes  are  placed  through  interpolation  if  the  spacing  between 
two  consecutive  control  volumes  is  greater  than  1  cm.  The 
collection  of  these  points  represents  the  ROI  contour  surface. 

The  contour  mapping  was  carried  out  in  three  steps:  (i) 
mapping  the  introduced  control  volumes  collectively  using  a 
rigid  image  registration  algorithm;  (ii)  iteratively  fine-tuning 
the  3D  positions  of  the  control  volumes  to  determine  the 
deformed  ROI  contour  surface;  and  (iii)  reconstructing  the 
new  contours  by  cutting  the  deformed  surface  slice-by-slice 
along  the  transversal,  sagittal,  or  coronal  direction.  This  pro¬ 
cedure  is  shown  in  Fig.  2  and  described  in  detail  below. 

II.D.  Collective  mapping  of  control  volumes 

After  a  series  of  control  volumes  were  placed  along  the 
segmented  contours  on  the  template  phase,  we  mapped  them 
onto  the  target  phase  collectively  (i.e.,  all  the  control  vol¬ 
umes  were  treated  as  an  entity)  using  a  rigid  image  registra¬ 
tion  algorithm.24  A  feature  of  the  rigid  collective  mapping  is 
that  the  relative  distances  and  orientations  of  the  control  vol¬ 
umes  remain  the  same  during  the  course  of  mapping,  result¬ 
ing  in  an  approximate  ROI  contour  surface  in  the  target 
phase  and  providing  a  good  start  for  further  adjustment  of 
the  contour  shape  to  accommodate  the  organ  deformations. 
We  note  that  using  a  rigid  mapping  is  not  a  necessary  step 
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Fig.  2.  Flow  chart  of  the  automatic  contour  mapping  process  for  4D  CT 
images. 


and  other  methods  capable  of  providing  a  reasonable  initial 
estimate  of  the  RO I  surface  may  also  be  used. 

The  collective  mapping  of  the  control  volumes  was  per¬ 
formed  under  the  guidance  of  a  normal  cross  correlation 
(NCC)  metric24  defined  by 

SXSUnW) 

f  =  — =  =  ,=1  ,  (i) 

V  a  M  p  j=l 

where  IJxj)  is  the  intensity  at  a  point  xj  in  a  control  volume 
indexed  by  a  in  the  target  phase  and  /^(x-)  represents  the 
intensity  at  a  point  x-  in  a  control  volume  indexed  by  /3  in 
the  template  phase.  In  Eq.  (1),  xand  x'  are  related  by  a  rigid 
transformation  T,  where  Tx'=x.  We  calculated  the  transfor¬ 
mation  matrix  T,  which  maps  the  ensemble  of  control  vol¬ 
umes  from  the  template  phase  to  the  target  phase.  The  lim¬ 
ited  memory  B  royden- FI  etcher- Goldfarb- Shannon 

algorithm  (L-BFGS)  algorithm26"30  was  used  to  optimize  the 
metric  in  Eq.  (1)  with  respect  to  the  transform  parameters. 
The  details  of  this  algorithm  have  been  described 
elsewhere18,31  and  will  not  be  repeated  here.  During  the 
course  of  the  control  volume  mapping,  an  iterative  calcula¬ 
tion  based  on  the  L-BFGS  algorithm  was  carried  out  until  a 
preset  maximum  number  of  iterations  were  met  or  until  the 
NCC  function  stopped  improving. 

II.E.  Fine  tuning  of  the  mapped  control  volumes 

In  the  absence  of  deformation,  the  RO  I  contour  of  the 
target  phase  is  obtained  by  connecting  the  centers  of  the 


rigidly  mapped  control  volumes.  In  a  more  general  case 
where  deformation  does  exist,  the  contour  from  the  rigid 
mapping  serves  as  an  initial  estimate  of  the  ROI  contour  in 
the  target  phase.  Further  positional  adjustment  of  the  control 
volumes  is  needed  to  accommodate  the  deformation  of  the 
ROI.  The  final  positions  of  the  control  volumes,  which  define 
the  ROI  surface,  are  determined  by  balancing  the  self- energy 
and  the  interaction  energy,  which  drive  the  control  volumes 
to  their  corresponding  locations  while  maintaining  similar 
shape  integrity  of  the  contour. 

Mathematically,  the  above  adaptation  process  is  modeled 
by  two  "energy  terms."  The  first  term  is  referred  to  as  "self¬ 
energy"  and  the  second  term  describes  the  "interaction" 
among  the  control  volumes.  Self-energy  tends  to  drive  a  con¬ 
trol  volume  toward  a  position  where  the  neighborhood  envi¬ 
ronment  resembles  itself  most.  For  each  control  volume  in 
the  template  image,  this  process  was  driven  by  the  NCC 
between  the  volume  and  its  corresponding  locations  in  the 
target  image.  The  interaction  term  among  the  control  vol¬ 
umes  intends  to  maintain  the  integrity  of  the  control  volume 
cluster  as  a  whole  and  prevents  any  unrealistic  control  vol¬ 
ume  configuration  from  happening.  The  interaction  energy 
term  is  described  by 

E  interaction  —  M  template  —  ^  target'  (2) 

where  M  template  is  defined  as  the  correlation  between  a  con¬ 
trol  volume  in  the  template  phase  and  the  neighboring  con¬ 
trol  volumes  in  the  target  phase,  M  target  represents  the  corre¬ 
lation  between  the  control  volume  in  the  target  phase  and  the 
neighboring  volumes  in  the  template  phase.  These  two  terms 
take  into  account  the  neighborhood  environment  of  the  con¬ 
trol  volume  being  adjusted  and  apply  a  constraint  on  the 
possible  form  of  the  control  volume  configuration.  In  a 
sense,  the  interaction  energy  in  Eq.  (2)  exerts  a  restoring 
force  when  the  position  of  a  control  volume  is  varied  with 
respect  to  its  neighbors.  The  interaction  force  is  important  in 
preventing  the  control  volumes  from  moving  to  unrealistic 
positions  simply  driven  by  the  self-energy  and  in  retaining 
the  shape  integrity  of  the  ROI  surface.  For  simplicity,  only 
the  adjacent  control  volumes  were  considered  when  comput¬ 
ing  the  interaction  energies.  The  final  position  of  each  con¬ 
trol  volume  was  determined  by  minimizing  the  sum  of  the 
self-energy  and  interaction  energy  terms. 

A  simple  searching  algorithm  was  implemented  to  find  the 
minimum  of  Eq.  (2)  and  thus  the  final  configuration  of  the 
control  volumes.  The  algorithm  was  realized  using  an  ex¬ 
haustive  search  in  the  defined  local  region  around  the  control 
volume  and  was  based  on  three  dimensional  basis.  For  each 
rigidly  mapped  control  volume  on  the  target  phase,  we  de¬ 
fined  a  small  region  of  ~3  cm  around  the  central  point  in 
the  volume.  The  search  region  of  3  cm  was  chosen  primarily 
to  accommodate  maximum  deformation  of  tissue.  3  cm  is  a 
very  conservative  value  because  the  tissue  deformation  in 
two  adjacent  phases  barely  exceeds  this  value  (the  maximum 
deformation  occurring  between  the  inhale  and  exhale  phases 
may  reach  ~3  cm).  Equation  (2)  was  then  minimized  to 
determine  the  optimal  control  volume  location  within  this 
region.  This  process  was  repeated  for  each  separate  control 
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volume.  Due  to  the  fact  that  the  interaction  energy  term  de¬ 
mands  the  information  from  the  neighboring  volumes,  ad¬ 
justing  the  location  of  each  control  volume  changes  the  po¬ 
sition  of  the  previously  fine-tuned  volume.  Therefore,  several 
cycles  of  this  process  were  needed  to  obtain  the  truly  optimal 
positions  of  all  the  control  volumes.  We  found  (see  digital 
phantom  study  in  results  section)  that  2-3  complete  adjust¬ 
ment  cycles  were  adequate  to  find  the  optimal  control  vol¬ 
ume  configuration. 

II.F.  Reconstruction  of  target  contours 

Upon  completion  of  the  control  volume  mapping  and  ad¬ 
aptation,  the  centers  of  each  control  volume  were  identified. 
A  ROI  surface  was  constructed  with  these  central  points  us¬ 
ing  a  triangulated  surface  construction  technique  which  uses 
marching  cubes32  method  and  triangular  surface  decimation 
of  VTK. 25,33  The  intersection  of  the  surface  with  each  CT 
slice  was  superimposed  on  top  of  the  image  in  order  to  vi¬ 
sualize  the  contour  in  a  conventional  fashion.  These  contours 
can  be  exported  as  ASCII  or  DICOM-RT  format  for  treat¬ 
ment  planning. 

II.G.  Evaluation  of  the  algorithm  and  case  study 

Evaluation  of  a  contour  mapping  algorithm  is  a  difficult 
task  because  of  the  general  lack  of  a  gold  standard  for  com¬ 
parison.  The  proposed  control  volume  based  4D  contouring 
technique  was  first  evaluated  with  a  digital  phantom  experi¬ 
ment.  In  this  study,  a  thoracic  CT  image  was  deformed  in¬ 
tentionally  using  a  known  deformation  matrix  which  was 
introduced  by  drifting  the  positions  of  each  control  volume 
along  the  contour  with  known  sizes.  Specifically,  for  a  con¬ 
trol  volume  /,  we  assign  the  following  displacements: 

D(x,y)  =/*  0.01  and  D(z)=i*  0.001  if  / 

=  0,  1 . N/2  (3) 

or 

D(x,y)  =(N  -/)*  0.01  and  D  (z)  =  (N  -/)*  0.001  if  /' 

=  A//2  +  1,  N/2  +2,  ...  A/,  (4) 

where  N  is  the  total  number  of  control  volumes  on  the  ROI 
contour.  The  artificially  deformed  image  serves  as  a  "new 
breathing  phase"  relative  to  the  original  image.  The  "ground 
truth"  lung  surface  in  the  target  phase  was  attainable  by 
transforming  the  original  contours  with  the  same  deforma¬ 
tion  matrix.  The  manually  delineated  lung  contours  in  the 
original  image  were  also  mapped  using  the  proposed  novel 
control  volume  based  approach.  A  comparison  of  the  mapped 
lung  surface  with  the  ground  truth  allowed  us  to  quantita¬ 
tively  assess  the  success  of  the  proposed  approach. 

The  control  volume  based  contour  mapping  technique 
was  also  applied  to  three  4D  CT  patient  scans.  A  physician 
manually  delineated  the  lungs  and  GTV  for  each  scan  on  a 
selected  phase  (for  example,  the  exhale  phase).  These  con¬ 
tours  were  then  mapped  onto  all  other  phases  using  the  pro¬ 
posed  technique.  Visual  inspection  was  used  to  evaluate  the 
three  patient  studies.  Although  it  is  less  quantitative,  visual 


(a)  (b) 


Fig.  3.  Validation  process  of  proposed  algorithm  using  a  digital  phantom, 
(a)  Original  image  and  contour;  (b)  Artificially  deformed  image  and  contour 
(ground  truth  contour)  together  with  contour  after  fine  tuning.  These  two 
contours  are  almost  indistinguishable.  The  contour  after  rigid  mapping  is 
also  shown. 


inspection  is  a  cnvenient  way  for  rapid  assessment  of  a  seg¬ 
mentation  calculation,  especially  in  a  case  where  the  ground 
truth  contours  do  not  exist. 

III.  RESULTS 

III.A.  Digital  phantom  study 

The  thoracic  CT  images  before  and  after  the  intentionally 
introduced  deformation  are  shown  in  Fig.  3.  The  manually 
delineated  lung  contour  is  shown  in  Fig.  3(a),  while  the 
ground  truth  contour  obtained  by  transforming  the  manual 
contour  with  the  known  transformation  matrix  is  plotted  in 
Fig.  3(b).  The  contour  from  the  proposed  approach  is  shown 
in  the  same  figure.  To  be  comprehensive,  the  rigidly  mapped 
contour,  which  served  as  the  start  of  the  model-based  refine¬ 
ment,  is  also  plotted  in  Fig.  3(b).  Clearly,  it  is  significantly 
deviated  from  the  lung  boundary.  Visual  inspection  of  the 
ground  truth  contour  and  contour  after  refinement  in  Fig. 
3(b)  indicated  that  the  gold  standard  and  the  mapped  con¬ 
tours  are  similar,  despite  the  fact  that  the  intentionally  intro¬ 
duced  deformation  field  is  quite  large  (the  displacement  of 
the  some  of  the  voxels  is  as  large  as  2.0  cm).  The  mean  and 
maximum  separations  between  the  two  sets  of  contours  were 
found  to  be  1.5  and  2.5  mm,  respectively. 

To  better  understand  the  mapping  process,  we  examined 
the  behavior  of  each  energy  term  (metrics)  during  the  course 
of  the  contour  mapping.  In  Fig.  4(a)  we  plotted  the  values  of 
self-energy,  interactive  energy,  and  total  energy  of  one  of  the 
control  volumes  in  the  process  of  positional  adjustment  in 
the  surrounding  area  of  the  control  volume  in  the  deformed 
image.  For  this  particular  control  volume,  the  minimum  of 
Eq.  (2)  was  reached  after  searching  the  first  30  points.  For 
control  volumes  located  in  regions  where  the  deformation  is 
large,  more  searching  points  may  be  required  to  reach  the 
minimum.  It  is  also  interesting  to  show  how  the  average 
metric  value  of  all  the  control  volumes  evolved.  In  Fig.  4(b) 
the  metric  value  is  depicted  as  a  function  of  the  calculation 
process.  The  first  point,  corresponding  to  step  0,  is  the  metric 
before  mapping.  Step  1  refers  to  the  system  after  rigid  map¬ 
ping.  The  metric  at  this  point  is  decreased  but  does  not  reach 
the  minimum.  The  third  point  shows  the  metric  after  each  of 
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Fig.  4.  (a)  Energies  (self,  interaction,  and  total  energy)  as  a  function  of  the  searching  points  in  the  surrounding  area  of  a  randomly  selected  control  volume 
in  the  target  image  (see  text  for  details),  (b)  The  metric  value  versus  step  during  the  course  of  the  contour  mapping.  The  metric  values  are  the  average  NCC 
of  all  the  involved  control  volumes. 


the  control  volumes  is  fine  tuned  sequentially.  To  obtain  the 
optimal  contours  two  more  cycles  of  fine  tuning  were  per¬ 
formed;  however,  the  improvement  with  further  iterations 
was  not  significant. 

III.B.  Patient  studies 

The  proposed  technique  was  assessed  using  4D  CT  im¬ 
ages  of  three  lung  cancer  patients.  Figure  5  shows  the  manu¬ 
ally  delineated  contours  on  the  4D  CT  scans  of  two  lung 
cancer  patients  for  lung  and  GTV,  respectively.  The  CT  im¬ 
ages  and  contours  are  displayed  in  axial,  coronal,  and  sagittal 
views.  Contours  were  drawn  on  the  exhale  phases  (phase  1) 
for  patient  1  (top  row  in  Fig.  5).  The  GTV  contour  for  patient 
3  was  delineated  on  the  inhale  phase  (phase  5)  as  shown  in 
Fig.  5  bottom  row.  For  patients  1  and  2  the  mapping  of  the 
lung  contours  was  studied.  For  the  third  patient,  GTV  con¬ 
tour  mapping  was  investigated. 

The  lung  contours  before  and  after  the  control  volume- 
based  mapping  for  the  first  patient  are  presented  in  Fig.  6.  To 


(a)  (b)  (c) 


Fig.  5.  The  template  CT  images  and  manually  segmented  contours  for  two 
lung  cancer  patients  (top  row:  patient  1;  bottom  row:  patient  3)  on  axial, 
coronal,  and  sagittal  views.  Right  lung  of  patient  1  and  GTV  of  patient  3 
were  manually  delineated,  respectively. 


better  visually  evaluate  the  results,  the  CT  images  and  con¬ 
tours  for  this  patient  are  displayed  in  axial,  coronal,  and  sag¬ 
ittal  views  as  in  left,  middle,  and  right  columns,  respectively. 
Target  phases  4,  7,  and  10  are  selected  to  demonstrate  the 
results  as  in  the  top,  middle,  and  bottom  rows  in  Fig.  6.  The 
template  contours  and  rigidly  mapped  contours  were  shown 
and  significantly  deviated  from  the  ROI  boundary.  The  final 
contours  after  the  model-based  adaption  are  presented  as 
well. 

When  the  lung  deformation  is  small,  for  example,  in 
phase  10  of  patient  1  [see  Figs.  6(g)- 6(i)],  the  rigidly 
mapped  contours  closely  resemble  the  target  contours.  Fine 
tuning  merely  provides  limited  adjustment  of  the  contours. 
For  phases  with  large  deformations  (e.g.,  phases  of  4  and  7 


Fig.  6.  Axial,  coronal,  and  sagittal  views  of  CT  images  along  with  lung 
contours  for  the  first  patient.  Top  row:  phase  4;  middle  row:  phase  7;  bottom 
row:  phase  10.  Rigidly  mapped  and  target  contours  are  displayed.  Template 
contour  is  also  overlaid  on  displayed  phases. 
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Fig.  7.  Axial,  coronal,  and  sagittal  views  of  CT  images  along  with  GTV 
contours  for  the  third  patient.  Top  row:  phase  2;  middle  row:  phase  4; 
bottom  row:  phase  7.  Rigidly  mapped  and  target  contours  are  displayed. 
Template  contour  is  also  overlaid  on  displayed  phases. 


of  this  patient),  the  model-based  adjustment  after  the  rigid 
contour  mapping  plays  an  important  role  in  obtaining  opti¬ 
mal  contour  arrangement.  In  case  of  the  lung  contours,  our 
analyses  of  the  patient  data  indicate  that  the  accuracy  of  our 
contour  mapping  technique  is  within  one  pixel  in  the  supe¬ 
rior  part  of  the  lung,  where  respiration  induced  deformation 
is  small.  In  other  parts,  the  average  error  on  the  target  con¬ 
tours  is  estimated  to  be  less  than  2.5  mm. 

The  mapped  GTV  contours  for  phases  2,  4,  and  7  in  the 
study  of  patient  3  are  shown  in  Fig.  7.  The  representation  of 
the  contours  is  the  same  as  in  Fig.  6.  For  phases  4  and  7  the 
deformation  was  relatively  small  and  our  contour  mapping 
algorithm  performed  reliably.  For  phase  2,  which  corre¬ 
sponded  to  the  exhale  phases,  significant  deformation  in  the 
ROIs  was  observed.  Despite  this  less  ideal  case,  our  algo¬ 
rithm  still  worked  well  in  these  phases. 

IV.  DISCUSSION 

As  in  3D  radiation  therapy,  delineation  of  ROIs  in  4D  CT 
images  is  a  necessary  step  for  treatment  planning.34  4D  seg¬ 
mentation  is  required  for  constructing  a  4D  patient  model 
and  for  computing  the  accumulated  dose  of  moving  or  de¬ 
formed  organs.  One  way  to  proceed  is  to  draw  all  contours 
on  one  of  the  phases  of  respiration  and  map  or  propagate 
these  contours  onto  the  remaining  phases.  Various  deform¬ 
able  models  have  been  used  for  ROI  contour  mapping  as 
discussed  in  the  introduction.  In  this  work,  we  take  a  re¬ 
gional  approach  based  on  the  mapping  and  adaptation  of  the 
sparsely  sampled  control  volumes.  Our  approach  takes  ad¬ 
vantage  of  the  imaging  features  surrounding  the  ROI  and 
uses  them  as  guidance  in  searching  for  the  optimal  mapped 
contours  while  considering  the  shape  integrity  of  the  ROI 
surface. 


The  mapping  of  a  point  in  one  image  to  another  is  easily 
achievable  if  a  unique  identifier  or  signature  can  be  tagged  to 
the  point.  FI  ere  the  image  feature  contained  in  a  control  vol¬ 
ume  is  employed  as  a  signature  of  the  point  to  facilitate  the 
process  of  finding  its  corresponding  location  in  the  target 
phases.  The  concept  of  control  volume  was  first  introduced 
by  Schreibmann  and  Xing23  and  its  advantages  for  both 
intra-  and  intermodality  image  registration  have  been  dem¬ 
onstrated.  This  study  represents  a  novel  application  of  the 
concept  to  4D  image  segmentation.  After  a  rigid  mapping  of 
the  contours  from  the  template  phase  to  the  target  phase,  a 
model-based  adaptation  is  performed  to  establish  a  reliable 
association  between  the  ROIs  in  two  phase  specific  image 
sets.  This  adaptation  takes  into  account  the  deformation  of  an 
object  and  ensures  the  overall  integrity  of  the  resultant  ROI 
surface.  The  calculation  is  local  in  nature,  which  improves 
both  computational  efficiency  and  convergence  behavior. 

A  quantitative  comparison  of  different  deformable  regis¬ 
tration  algorithms  is  a  difficult  task  because  of  the  multifac¬ 
eted  and  even  subjective  nature  of  the  problem.  A  unbiased 
and  meaningful  comparison  may  entail  the  efforts  from  mul¬ 
tiple  institutions.35  Thus  we  defer  the  detailed  comparative 
study  to  the  future.  Flowever,  we  wish  to  emphasize  that  the 
chief  advantage  of  the  proposed  technique  is  its  computa¬ 
tional  efficiency.  Because  of  the  regional  nature  of  the  calcu¬ 
lation,  our  general  finding  is  that  the  new  algorithm  is  at  least 
an  order  of  magnitude  faster  than  the  whole  image  based 
approach.17,18  Enormous  saving  in  the  memory  usage  in  the 
proposed  approach  is  also  self-explanatory  and  highly  desir¬ 
able  feature  in  practice. 

A  common  problem  in  image  segmentation  and  contour 
mapping  studies  is  the  lack  of  quantitative  validation.  In  the 
study  of  Lu  et  a  I.,36  for  example,  the  accuracy  of  a  deform¬ 
able  model-based  contour  mapping  technique  was  evaluated 
purely  based  on  visual  inspection.  The  same  approach  was 
employed  in  many  other  previous  investigations.4,6,16,18  In 
our  study,  in  addition  to  the  visual  evaluation,  a  set  of  digital 
phantom  experiments  was  introduced  to  evaluate  the  success 
of  the  proposed  technique.  By  applying  a  prespecified  defor¬ 
mation  matrix  to  the  original  image,  the  ground  truth  of  the 
contour  propagation  is  readily  known.  Therefore,  the  experi¬ 
ments  provide  a  quantitative  test  of  the  proposed  algorithm. 
In  general,  we  found  that  a  spatial  accuracy  better  than  2.5 
mm  is  achievable  using  our  technique. 

V.  CONCLUSION 

The  development  of  4D  radiation  therapy  involves  the  use 
of  a  large  number  of  images  acquired  at  different  times 
and/or  with  different  modalities.  Clinical  implementation  of 
the  new  IGRT  paradigm  is,  to  a  large  extent,  bottlenecked  by 
the  inability  to  accurately  and  efficiently  register  images  and 
segment  the  ROIs.  In  this  work,  a  mathematical  framework 
for  control  volume-based  contour  mapping  has  been  pro¬ 
posed  for  4D  radiation  therapy.  We  demonstrated  that  the 
information  contained  in  the  boundary  region  is  sufficient  to 
guide  the  contour  mapping  process  without  registering  the 
whole  image  or  relying  on  the  use  of  an  ad  hoc  surface 
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deforming  model.  The  results  showed  that  the  control 
volume-based  contour  mapping  algorithm  is  capable  of  ro¬ 
bustly  and  accurately  mapping  contours  from  one  phase  of  a 
4D  CT  to  the  remaining  phases.  Our  technique  decreases  the 
workload  involved  in  4D  CT  ROI  segmentation  and  provides 
a  valuable  tool  for  the  efficient  use  of  available  spatial -tempo 
information  for  4D  simulation  and  treatment. 
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AUTOMATED  CONTOUR  MAPPING  WITH  A  REGIONAL  DEFORMABLE  MODEL 
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Albert  Koong,  M.D.,  and  Lei  Xing,  Ph.D. 

Department  of  Radiation  Oncology,  Stanford  University  School  of  Medicine,  Stanford,  CA 

Purpose:  To  develop  a  regional  narrow-band  algorithm  to  auto-propagate  the  contour  surface  of  a  region  of  inter¬ 
est  (ROI)  from  one  phase  to  other  phases  of  four-dimensional  computed  tomography  (4D-CT). 

Methods  and  Materials:  The  ROI  contours  were  manually  delineated  on  a  selected  phase  of  4D-CT.  A  narrow  band 
encompassing  the  ROI  boundary  was  created  on  the  image  and  used  as  a  compact  representation  of  the  ROI  sur¬ 
face.  A  BSpline  deformable  registration  was  performed  to  map  the  band  to  other  phases.  A  Mattes  mutual  infor¬ 
mation  was  used  as  the  metric  function,  and  the  limited  memory  Broyden-Fletcher-Goldfarb-Shanno  algorithm 
was  used  to  optimize  the  function.  After  registration  the  deformation  field  was  extracted  and  used  to  transform 
the  manual  contours  to  other  phases.  Bidirectional  contour  mapping  was  introduced  to  evaluate  the  proposed  tech¬ 
nique.  The  new  algorithm  was  tested  on  synthetic  images  and  applied  to  4D-CT  images  of  4  thoracic  patients  and 
a  head-and-neck  Cone-beam  CT  case. 

Results:  Application  of  the  algorithm  to  synthetic  images  and  Cone-beam  CT  images  indicates  that  an  accuracy  of 
1.0  mm  is  achievable  and  that  4D-CT  images  show  a  spatial  accuracy  better  than  1.5  mm  for  ROI  mappings  be¬ 
tween  adjacent  phases,  and  3  mm  in  opposite-phase  mapping.  Compared  with  whole  image-based  calculations, 
the  computation  was  an  order  of  magnitude  more  efficient,  in  addition  to  the  much-reduced  computer  memory 
consumption. 

Conclusions:  A  narrow-band  model  is  an  efficient  way  for  contour  mapping  and  should  find  widespread  applica¬ 
tion  in  future  4D  treatment  planning.  ©  2008  Elsevier  Inc. 

Deformable  model,  Image  registration,  Contour  mapping,  IGRT. 


INTRODUCTION 

Segmentation  of  a  region  of  interest  (ROI),  such  as  a  tumor 
target  volume  or  a  sensitive  structure,  is  an  important  but 
time-consuming  task  in  radiotherapy  (1-6).  With  the  emer¬ 
gence  of  four-dimensional  (4D)  imaging  and  adaptive  radio¬ 
therapy,  the  need  for  efficient  and  robust  segmentation  tools 
is  even  increasing  (7-11).  Because  of  dramatically  increased 
numbers  of  images,  it  becomes  impractical  to  manually  seg¬ 
ment  the  ROIs  slice  by  slice  as  in  current  three-dimensional 
radiotherapy  practice.  A  natural  solution  to  the  4D  computed 
tomography  (4D-CT)  segmentation  problem  is  to  delineate 
the  ROIs  on  a  selected  phase  and  then  propagate  the  contours 
onto  other  phases  using  a  mathematical  model.  Along  this 
line,  deformable  model-based  contour  mapping  has  been 
implemented  by  a  few  groups  (12-14).  Although  feasible, 
the  calculation  is  global  in  nature  and  thus  computationally 
intensive.  In  addition,  the  accuracy  of  the  mapped  contours 
may  be  compromised  because  the  registration  may  be 


influenced  unnecessarily  by  the  image  content  distant  from 
the  ROIs,  which  would  otherwise  be  irrelevant  to  the  contour 
mapping  process.  This  is  especially  problematic  when  non¬ 
local  deformable  models,  such  as  thin  plate  spline  and  elastic 
model,  are  used.  In  general,  contour  mapping  is  a  regional 
problem,  and  a  global  association  of  the  phase-based  images 
is  neither  necessary  nor  efficient. 

Surface  mapping  techniques  (15-17)  represent  a  competi¬ 
tive  alternative  to  the  deformable  model-based  approach. 
The  idea  of  surface  mapping  is  to  obtain  contour  transforma¬ 
tion  by  iteratively  deforming  the  ROI  contour-extended  sur¬ 
face  until  the  optimal  match  with  the  reference  is  found.  The 
calculation  involves  only  the  surface  region  and  is  thus  com¬ 
putationally  efficient.  Numerous  surface  mapping  techniques 
have  been  developed  in  the  past,  which  include,  to  name 
a  few,  spatial  partitioning,  principal  component  analysis, 
conformal  mapping,  rigid  affine  transformation,  deformable 
contours,  and  warping  based  on  the  thin-plate  spline.  All  of 
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these  techniques  are  a  mapping  between  topologic  compo¬ 
nents  of  the  input  surfaces  that  allow  for  transfer  of  annota¬ 
tions.  Although  the  calculations  are  inherently  efficient,  the 
results  depend  heavily  on  the  model  used,  which  may  not 
be  generally  applicable  for  all  clinical  situations  because 
the  ROI  surface  is  multidimensional  and  hardly  modeled 
by  only  a  few  parameters. 

In  this  work,  we  present  a  novel  regional  algorithm  for  ROI 
propagation  among  different  4D-CT  phases.  The  deforma¬ 
tion  of  an  ROI  contour-extended  surface  in  our  algorithm  is 
not  driven  by  an  ad  hoc  surface-based  model  but  instead  by 
the  image  features  in  the  neighborhood  of  the  surface.  The 
underlying  hypothesis  here  is  that  information  contained  in 
the  ROI  boundary  region  is  sufficient  to  guide  the  contour 
mapping  process.  In  the  proposed  algorithm  the  neighbor¬ 
hood  image  features  of  an  ROI  are  captured  by  a  narrow 
band,  which  is  composed  of  all  points  within  two  surfaces 
with  the  signed  distances  of  =bd  from  the  ROI  boundary. 
The  algorithm  is  a  hybrid  of  the  regional  surface-based 
model  and  the  global  deformable  registration-based  ap¬ 
proach.  The  combination  takes  advantage  of  the  desirable 
features  of  each  of  these  two  techniques  and  provides  a  robust 
and  computationally  efficient  contour  propagation  tool  for 
4D  radiotherapy. 

METHODS  AND  MATERIALS 

Software  platform 

The  proposed  contour  mapping  algorithm  was  implemented  using 
the  Insight  Toolkit  (18)  and  the  Visualization  Toolkit  (19),  which 
are  open  source  cross-platform  C++  software  toolkits  sponsored 
by  the  National  Library  of  Medicine. 

Overview  of  the  mapping  process 

Figure  1  depicts  the  overall  contour  mapping  process.  For  a  given 
4D-CT  image  set,  a  selected  phase,  named  the  template  phase,  was 
selected,  and  the  ROIs  were  manually  delineated  by  a  physician.  The 
manually  outlined  contour  was  referred  to  as  the  template  contour.  A 
narrow  band  encompassing  the  template  contour  was  created  (see 
next  section  for  details).  A  deformable  mapping  was  then  carried 
out  to  propagate  the  band  from  the  template  phase  to  other  phases, 
referred  to  as  target  phases.  Upon  successful  mapping  of  the  band, 
the  deformation  field  was  used  to  transform  the  template  contour 
to  the  target  images. 

Narrow -band  representation  of  ROI  contour 

The  contour  manually  segmented  on  an  axial  slice  of  the  template 
image  has  a  polygon  shape,  and  the  vertices  of  the  polygon  form  the 
basis  for  constructing  the  narrow  band.  As  schematically  shown 
in  Fig.  2,  a  band  with  signed  distances  ±d  was  placed  along  the  tem¬ 
plate  contour.  The  regional  image  features  contained  in  the  band 
function  serve  as  a  “signature”  of  the  contour  and  drive  the  contour 
mapping  process.  The  distance  between  the  neighboring  vertices  on 
the  contour  is  typically  2-10  mm,  depending  on  the  shape  of  the 
contour.  In  generating  the  narrow  band,  we  first  created  cubes 
with  a  side  length  of  2d  around  all  the  vertices,  as  depicted  by  points 
A  and  B  in  Fig.  2.  To  obtain  a  smooth  band,  between  A  and  B  three 
more  cubes,  centered  at  points  C,  D,  and  E,  were  inserted.  Point  C 
was  chosen  to  be  the  middle  point  between  A  and  B,  point  D  the 
middle  between  A  and  C,  and  point  E  the  middle  between  C  and 
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(a) 


(b) 


Fig.  1.  Flow  chart  of  narrow  band-based  contour  mapping  proce¬ 
dure.  (a)  Overall  calculation  process,  (b)  Deformable  mapping  pro¬ 
cess  of  the  narrow  band. 

B.  More  interpolated  vertex  points  can  be  introduced  similarly 
when  needed.  Figure  3  illustrates  a  narrow  band  surrounding  the 
lung  boundary  on  the  template  phase  CT  image.  The  light  green 
area  stands  for  the  narrow  band,  and  the  green  curve  is  the  manual 
contour.  The  width  of  the  narrow  band  was  set  to  be  2d  =  15  mm 
in  our  calculations.  To  examine  the  robustness  of  the  proposed  map¬ 
ping  algorithm,  a  variety  of  other  bandwidths,  ranging  from  4  mm 
through  30  mm,  were  also  tested  for  one  of  the  clinical  cases. 

Contour  propagation 

As  illustrated  in  Fig.  1,  the  process  of  contour  mapping  is  essen¬ 
tially  to  warp  the  narrow  band  constructed  above  in  such  a  way 
that  its  best  match  in  the  target  image  is  found.  Mathematically, 
the  mapping  process  of  the  narrow  band  constitutes  an  optimization 
problem,  in  which  a  group  of  transformation  parameters  that  trans¬ 
form  the  points  within  the  band  in  the  template  phase  to  their  homol¬ 
ogous  points  in  the  target  image.  The  warping  of  the  narrow  band  is 
quantified  by  a  metric  function,  which  ranks  a  trial  matching  based  on 
the  “accordance”  level  of  the  image  content  of  the  band  and  its  cor¬ 
respondence  in  the  target  image.  The  calculation  process  is  detailed 
below. 
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Narrow  Band 


Fig.  2.  A  schematic  drawing  of  narrow-band  construction. 


The  input  to  the  contour  mapping  software  includes  the  narrow 
band  and  the  whole  target  image,  which  are  described  by  the  image 
intensity  distributions  Ia(x)  and  Ib(x),  respectively.  It  is  worth  em¬ 
phasizing  that,  even  though  the  whole  target  image  was  used,  only 
fractional  voxels  in  the  target  image  (the  voxels  encompassed  by 
the  band)  are  involved  in  each  iteration  (a  subregion  surrounding 
the  ROI  on  the  target  image  could  be  created  and  used  in  the  calcu¬ 
lation,  but  the  algorithm  converged  so  fast  that  after  two  to  three  it¬ 
erations  the  searching  was  quickly  confined  in  the  neighborhood  of 
the  optimal  solution).  The  narrow  band  acts  as  a  representation  of 
the  ROI  contour.  The  task  is  to  find  the  transformation  matrix, 
r(x),  that  maps  an  arbitrary  point  in  the  band  to  the  corresponding 
point  on  the  target  image  (or  vice  versa)  so  that  the  best  possible  cor¬ 
respondence,  as  measured  by  the  metric  function,  is  achieved.  The 
calculation  proceeds  iteratively.  A  B  Spline  deformable  model  is 
used  to  model  the  deformation  of  the  band,  but  other  models  should 
also  be  applicable.  The  spacing  between  the  BSpline  nodes  was  cho¬ 
sen  to  be  approximately  0.5  cm  (smaller  spacing  was  tested,  but  no 
significant  difference  was  found  in  the  final  registration  results).  The 
displacement  of  a  node  i  is  specified  by  a  vector  xi?  and  the  displace¬ 
ment  vectors  (20)  of  a  collection  of  nodes  characterize  the  tissue 
deformation.  The  displacement  at  a  location  x  on  the  image  is  de¬ 
duced  by  a  BSpline  polynomial  fitting. 

The  Mattes  Mutual  Information  (MMI)  (21)  was  used  as  the  met¬ 
ric  function  for  narrow-band  mapping  (22-25).  The  central  concept 


(C) 


Fig.  3.  Computed  tomographic  images  with  manual  contours  and  the  narrow  bands  for  patient  1.  The  narrow  bands  are 
shown  in  light  green  and  the  contours  are  green  curves,  (a)  Transverse  view;  (b)  coronal  view;  (c)  sagittal  view. 
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of  mutual  information  (MI)  is  the  calculation  of  entropy.  For  an  im¬ 
age  A,  the  entropy  is  defined  as 


H  A 


=  -  J pA(a)log  pA(a)da, 


where  p/fa)  (also  called  the  marginal  probability  density  function 
[PDF])  is  the  probability  distribution  of  grey  values  (image  intensi¬ 
ties),  which  is  estimated  by  counting  the  number  of  times  each  grey 
value  occurs  in  the  image  and  dividing  those  numbers  by  the  total 
number  of  occurrences.  Given  two  images,  A  and  B,  their  joint  en¬ 
tropy  is 


=  -  I  I  pAB(a,  b) log  pAB(a,  b)dadb, 

where  PAB(a,b)  is  the  joint  PDF  defined  by  a  ratio  between  the 
number  of  grey  values  in  the  joint  histogram  (feature  space)  of 
two  images  and  the  total  entries  (26).  The  mutual  information  is  gen¬ 
erally  expressed  as 

MI(A,B)  =  H(A)  +H(B)  -H(A,B). 

Mutual  information  measures  the  level  of  information  that  a  ran¬ 
dom  variable  (e.g.,  Ia(x))  can  predict  about  another  random  variable 
(e.g.,  Ib(x)).  Different  from  the  conventional  MI,  whereby  two  sep¬ 
arate  intensity  samples  are  drawn  from  the  image,  the  Mattes  imple¬ 
mentation,  MMI,  uses  only  one  set  of  intensity  to  evaluate  both  the 
marginal  and  joint  PDFs  at  discrete  positions  or  bins  that  uniformly 
spread  within  the  dynamic  range  of  the  images.  Entropy  values  were 
computed  by  summing  over  all  the  bins.  The  number  of  bins  used  to 
compute  the  entropy  in  MMI  metric  evaluation  was  chosen  to  be  30, 
and  the  number  of  spatial  samples  used  was  20,000.  Details  of  MMI 
implementation  can  be  found  in  Mattes  et  al.  (21). 

The  limited  memory  Broyden-Fletcher-Goldfarb-Shanno  algo¬ 
rithm  (L-BFGS)  (27-29)  was  used  to  optimize  the  MMI  metric  func¬ 
tion  with  respect  to  the  displacement  parameters  of  the  nodes,  { xA} , 
to  find  the  transformation  matrix  T(x)  that  relates  the  points  on 
image  A  and  image  B.  Here  we  just  briefly  show  the  algorithm. 
Starting  from  a  positive  definitive  approximation  of  the  inverse 
Hessian  H0  at  x0,  L-BFGS  derives  the  optimization  variables  by 
iteratively  searching  through  the  solution  space.  At  an  iteration  k, 
the  calculation  proceeds  as  follows:  [1]  determine  the  descent 
direction  =  —  HkVf(xk);  [2]  line  search  with  a  step  size 
otk  =  (xfc  +  op*),  where  a  is  the  step  size  defined  in  the  L- 

BFGS  software  package;  [3]  update  xk+1  =  xk  +  ak  pk ;  and  [4]  com¬ 
pute  Hk+1  with  the  updated  Hk  . 

At  each  iteration  a  backtracking  line  search  is  used  in  L-BFGS  to 
determine  the  step  size  of  movement  to  reach  the  minimum  of  / 
along  the  ray  xk  +  a:pk.  For  convergence  a  has  to  be  chosen  such 
that  a  sufficient  decrease  criterion  is  satisfied,  which  depends  on 
the  local  gradient  and  function  value  and  is  specified  in  L-BFGS 
by  the  Wolfe  conditions  (27).  During  the  course  of  optimization, 
the  above  iterative  calculation  based  on  L-BFGS  algorithm  con¬ 
tinues  until  the  following  stopping  criterion  is  fulfilled: 

jlfyfr*)  n2  <£ 

max(l,  \\xk\\2) 

or  a  pre-set  maximum  number  of  iterations  is  reached.  In  this  study 
we  set  e  =  106  and  the  iteration  number  to  200,  but  no  more  than  100 
iterations  were  exceeded  in  all  our  calculations  for  the  algorithm  to 
converge. 


Evaluation  of  algorithm  performance 

Evaluation  of  a  contour  mapping  algorithm  is  a  difficult  task  be¬ 
cause  of  the  lack  of  the  ground  truth  for  comparison.  A  straightfor¬ 
ward  means  of  evaluation  is  the  visual  inspection  of  the  mapped 
contours.  In  addition  to  this,  evaluation  based  on  synthetic  images 
(digital  phantoms)  is  also  commonly  used.  The  images  and  existing 
contours  are  distorted  with  preset  deformation  fields.  Because  the 
gold  standard  is  known,  a  direct  comparison  with  the  mapped  con¬ 
tour  is  made  so  as  to  assess  the  propagation  algorithm  quantitatively. 
Beside  these  two  methods,  we  further  performed  a  bidirectional 
mapping  to  evaluate  the  proposed  algorithm.  In  this  test,  the  reverse 
of  the  original  contour  mapping  was  performed:  the  mapped  con¬ 
tours  on  the  target  phase  were  treated  as  the  template  contours  and 
mapped  back  to  the  original  template  phase.  The  contours  so  ob¬ 
tained  were  then  compared  with  the  original  manual  contours,  and 
the  difference  between  the  two  sets  of  contours  was  quantified. 
The  difference  between  the  resultant  and  template  contours  was 
measured  in  terms  of  the  displacements  of  the  vertex  points  on  the 
two  contours.  The  last  yet  pragmatic  evaluation  of  the  algorithm  per¬ 
formance  on  patient’s  study  was  based  on  the  physician’s  manual 
contours. 

Case  study 

Four  thoracic  cancer  patients,  named  as  patient  1 , 2, 3,  and  4,  were 
first  used  to  test  the  proposed  algorithm.  These  patients  underwent 
4D-CT  scans.  The  4D-CT  images  were  acquired  with  a  GE  Discov- 
ery-ST  CT  scanner  (GE  Medical  System,  Milwaukee,  WI).  The  col¬ 
lected  data  were  sorted  into  10  phase  bins.  The  ROIs  on  the  template 
phase  were  manually  segmented  by  a  physician.  Specifically,  for  pa¬ 
tients  1  and  2,  the  inhale  phase  was  chosen  for  manual  segmentation, 
and  for  patients  3  and  4,  the  exhale  phase.  Different  ROIs  were  used 
to  better  evaluate  the  algorithm.  Lungs  were  selected  from  patients 
1,  2,  and  3  and  gross  tumor  volume  (GTV)  from  patient  4.  Figure  3 
illustrates  the  manual  contour  and  narrow  band  representation  for 
the  lung  from  patient  1 .  Contour  is  shown  in  the  green  curve  together 
with  the  regional  narrow  bands  (light  green  area)  on  the  transverse, 
coronal,  and  sagittal  views  (Figs.  3a,  3b,  and  3c,  respectively). 

To  further  assess  the  robustness  of  the  proposed  algorithm,  we 
also  carried  out  the  contour  propagation  calculation  from  planning 
CT  to  Cone-beam  CT  (CBCT)  for  a  head-and-neck  case.  The 
CBCT  images  were  acquired  using  the  Varian  Trilogy  system  (Var- 
ian  Medical  Systems,  Palo  Alto,  CA). 

RESULTS 


Convergence  analysis 

To  better  illustrate  the  iterative  process  of  the  contour 
propagation,  in  Fig.  4  the  MMI  metric  as  a  function  of  itera¬ 
tion  step  is  plotted  for  the  narrow  band  mapping  from  the  first 
phase  (inhale  phase)  to  the  other  nine  phases  for  the  first  tho¬ 
racic  patient.  In  all  nine  calculations  it  is  seen  that  the  metric 
value  decreases  monotonically  as  the  iteration  proceeds. 
However,  the  number  of  iterations  needed  for  the  algorithm 
to  find  the  optimal  solution  varies.  It  is  interesting  to  observe 
that,  for  an  “easier”  mapping  whereby  the  deformation  be¬ 
tween  the  two  phases  is  small,  the  number  of  iterations 
required  is  less,  whereas  for  “tougher”  ones  with  larger  dif¬ 
ferences  in  ROI  shapes,  the  required  number  of  iterations  in¬ 
creases  drastically.  Indeed,  from  Fig.  4  it  is  seen  that  the 
minimum  number  of  iterations  required  for  the  metric  to  sat¬ 
urate  occurs  when  mapping  the  phase  1  to  the  adjacent 
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Iteration 

Fig.  4.  Narrow-band  metric  values  as  a  function  of  iteration  step 
when  mapping  the  narrow  band  from  phase  1  to  the  other  nine 
phases  of  the  four-dimensional  computed  tomography. 

phases,  2  and  10.  For  other  mapping,  the  required  iteration  in¬ 
creases  and  reaches  its  largest  value  for  the  “toughest”  map¬ 
ping  between  inhale  and  exhale  (phase  5)  phases. 

In  the  above  analysis,  the  bandwidth  was  set  to  be  15  mm. 
The  performance  of  the  proposed  algorithm  was  also  evalu¬ 
ated  by  varying  the  width  in  the  range  of  4  mm  and  30 
mm.  Specifically,  we  tried  the  widths  of  4  mm,  8  mm,  10 
mm,  15  mm,  20  mm,  and  30  mm.  Our  results  revealed  that, 
when  the  band  was  too  narrow  (e.g.,  4  mm),  the  mapping 
may  fail  locally  at  a  place  not  containing  sufficient  neighbor¬ 
hood  image  features.  The  situation  is  improved  dramatically 
as  the  bandwidth  increases.  For  all  the  clinical  cases  studied 
here,  no  single  failure  was  observed  for  a  width  of  15  mm. 
When  the  width  is  too  large,  the  whole  ROI  will  be  included 
in  the  band.  In  this  situation,  the  mapping  becomes  equiva¬ 
lent  to  registering  the  whole  image  and  the  advantage  of 
the  narrow  band  will  be  overshadowed  by  the  dramatically 
increased  memory  and  computing  costs.  Our  experience  indi¬ 
cates  that  a  width  of  10-15  mm  provides  a  fine  balance 
between  the  computational  accuracy  and  the  associated  cost. 

We  found  that  the  overall  computing  time  was  increased 
by  roughly  an  order  of  magnitude  when  going  from  the  nar¬ 
row  band  approach  to  the  conventional  deformable  model- 
based  contour  mapping,  say,  approximately  3  min  for  narrow 
band-based  mapping  vs.  approximately  25  min  for  whole 
image-based  mapping.  The  dramatically  increased  computer 
memory  requirement  in  the  latter  case  also  posts  a  serious 
problem  when  developing  a  clinically  practical  contour  prop¬ 
agation  method  for  4D  radiotherapy. 

Algorithm  performance  evaluation 

In  addition  to  visual  inspect,  the  proposed  algorithm  was 
assessed  by  a  series  of  synthetic  images  or  digital  phantoms. 
Typically,  a  thoracic  CT  image  together  with  the  contour  was 
distorted  with  the  intentionally  introduced  deformation,  and 
then  the  contour  was  propagated  onto  the  distorted  image. 


A  quantitative  comparison  was  carried  out.  The  mean  and 
maximum  separation  between  the  gold  standard  and  the  map¬ 
ped  contours  were  found  to  be  1.0  mm  and  1.5  mm,  respec¬ 
tively.  Figure  5  shows  one  example  of  digital  phantom 
experiments. 

The  performance  of  the  proposed  algorithm  was  further 
evaluated  by  the  bidirectional  mapping  calculation  outlined 
in  Methods  and  Materials.  A  template  contour  at  phase  1 
was  first  mapped  to  phases  3  and  6.  The  mapped  contours 
were  then  treated  as  the  “starting  contours”  and  mapped 
back  to  phase  1 .  The  two  back-mapped  contours  were  com¬ 
pared  with  the  original  template  contour.  The  displacement 
of  each  back-mapped  vertex  point  relative  to  its  original  loca¬ 
tions  was  computed,  and  a  mean  value  of  0.8  mm  was  found 
for  the  bidirectional  mapping  between  phases  1  and  3  and  1.8 
mm  between  phases  1  and  6.  The  larger  displacement  in  the 
latter  situation  was  due  to  the  fact  that,  computationally,  it  is 
more  difficult  to  map  between  two  opposite  phases,  such  as 
inhale  and  exhale  phases,  owing  to  larger  organ  deforma¬ 
tions.  Overall,  the  observed  displacement  is  comparable  to 
the  pixel  size,  indicating  that  the  mapping  is  accurate  and 
robust. 

Thoracic  patient  study  results 

Figure  6  shows  the  contour  mapping  results  for  the  first 
clinical  case.  The  results  are  presented  in  axial,  coronal, 
and  sagittal  planes  for  phases  2  (Fig.  6a-c),  6  (Fig.  6d-f),  8 
(Fig.  6g-i),  and  10  (Fig.  6j-l).  For  phases  2  and  10,  which 
are  immediately  adjacent  to  the  inhale  phase,  the  deformation 
is  relatively  small  and  the  mapped  contours  conform  to  the 
ROI  boundary  very  well.  This  represents  the  “easy”  map¬ 
ping  situation  and  is  consistent  with  the  analysis  presented 
above.  The  average  error  was  less  than  1.5  mm.  For  a  “re¬ 
mote”  phase,  such  as  phase  6  shown  in  Fig.  6d-f,  more 


Fig.  5.  Synthetic  image  and  overlaid  contours.  The  original  contour 
is  depicted  in  green,  gold  standard  contour  in  blue,  and  the  mapped 
contour  in  red. 
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Fig.  6.  Computed  tomographic  images  and  mapped  contours  for  thoracic  patient  1.  Displayed  are  selected  phases.  From 
the  top  row  to  bottom,  phases  2,  6,  8,  and  10  are  presented,  respectively.  For  each  phase,  transverse,  coronal,  and  sagittal 
views  are  shown  from  left  to  right. 


iterations  were  entailed  to  find  the  optimal  solution,  and  the 
resultant  contours  tend  to  be  worse  as  compared  with  those 
phases  adjacent  to  phase  1.  According  to  the  bidirectional 
mapping,  the  average  mapping  error  for  phase  6  was  esti¬ 
mated  to  be  less  than  3  mm.  The  mapped  GTV  contours  (in 
red)  together  with  manual  contours  (in  blue)  by  a  physician 


for  phases  1,  4,  8,  and  10  in  the  study  of  patient  4  are  shown 
in  Fig.  7  (parts  a,  b,  d,  and  e,  respectively).  The  template 
phase  (phase  6)  with  the  template  manual  contour 
(in  green)  is  shown  in  Fig.  7c.  In  addition,  the  template  man¬ 
ual  contour  from  this  phase  was  overlaid  on  all  the  displayed 
phases.  For  phases  4  and  8  the  deformation  was  relatively 
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Fig.  7.  Axial  view  of  computed  tomographic  images  with  gross  tumor  volume  contours  for  the  fourth  thoracic  patient,  (a), 
(b),  (c),  (d),  and  (e)  correspond  to  phases  1, 4, 6  (template  phase),  8,  and  10,  respectively.  The  green  curves  are  the  manually 
outlined  template  contour  from  phase  6,  and  the  red  curves  represent  the  contours  after  warping.  The  manual  contours  (in 
blue)  by  a  physician  on  individual  phases  were  also  displayed. 


small  (the  manual  contour  was  delineated  on  phase  6),  and 
fewer  iterations  were  needed  to  find  the  optimal  bands  on 
the  target  images.  For  phases  1  and  10,  whereby  deformation 
was  significant  in  the  ROIs  although  more  computing  load 
was  necessary,  a  good  result  was  still  achieved  with  our  nar¬ 
row-band  technique.  Comparisons  between  the  mapped  con¬ 
tours  and  the  manually  segmented  contours  by  physicians  for 
these  patients  were  also  performed,  and  results  revealed 
a  similar  level  of  accuracy  (maximum  and  mean  values  of 
the  discrepancy  between  the  two  sets  of  contours  are  2.8 
mm  and  -0.9  mm,  respectively). 

As  a  useful  application  of  the  proposed  technique,  in  Fig.  8 
we  present  the  mean  and  maximum  lung  displacements  of 
contour  vortices  for  each  breathing  phase  relative  to  their  lo¬ 
cations  on  the  template  phase.  As  seen  in  Fig.  8,  the  overall 
behavior  of  the  mean  and  maximum  displacements  is  consis¬ 
tent  with  our  intuitive  expectation.  For  cases  1  and  2,  the  in¬ 
hale  phase  (phase  1)  was  manually  segmented,  thus  the 
displacement  for  that  phase  is  zero.  For  other  phases,  both 
mean  and  maximum  displacement  values  vary  with  the 
breathing  phase  and  reach  their  maxima  at  the  opposite 
phase.  For  case  3  the  exhale  phase  was  manually  segmented, 
and  the  behavior  was  thus  opposite  to  cases  1  and  2.  In  gen¬ 
eral,  an  average  displacement  of  approximately  3  mm  was 
found  for  inhale  and  exhale  phases.  A  slight  digression  is  no¬ 
ticed  in  phase  7  of  patient  1,  which  may  be  caused  by  4D-CT 
binning  artifacts.  This  type  of  data  is  particularly  useful  in 
determining  the  patient-specific  tumor  margin  to  account 
for  breathing  motion  of  the  tumor  target. 


Contour  propagation  in  a  head-and-neck  case 

The  results  of  contour  mapping  for  the  head-and-neck  case 
are  summarized  in  Fig.  9.  Figure  9a  shows  the  planning  CT 
along  with  manually  delineated  contours,  and  Fig.  9b  dis¬ 
plays  the  mapped  contours  of  the  body,  mandible,  and 
GTV  on  CBCT.  For  body  and  mandible  a  simple  rigid  map¬ 
ping  is  enough  to  achieve  high  accuracy.  For  the  GTV,  how¬ 
ever,  the  proposed  deformable  registration  model  was 
necessary  to  adequately  propagate  the  contour.  A  visual 
inspection  of  the  propagated  contours  suggests  that  the  map¬ 
ping  is  clinically  acceptable. 

DISCUSSION 

Four-dimensional  CT  image  segmentation  represents 
a  necessary  step  in  constructing  a  4D  patient  model  and  com¬ 
puting  the  accumulated  dose  in  4D  radiotherapy.  A  natural 
way  to  tackle  the  problem  is  to  auto-map  the  manually  delin¬ 
eated  contours  on  one  of  the  phases  to  the  remaining  phases. 
In  this  work,  a  regional  computing  algorithm  was  introduced 
to  deal  with  the  issue.  The  approach  relies  on  the  assumption 
that  a  narrow  band  surrounding  the  manually  segmented  con¬ 
tour  can  capture  sufficient  information  to  drive  the  finding  of 
its  counterparts  in  other  phases  of  the  4D-CT.  Obviously,  this 
assumption  is  valid  when  the  band  is  sufficiently  wide  so  that 
a  large  number  of  voxels  are  involved  in  the  registration  cal¬ 
culation.  As  demonstrated  by  the  presented  data,  the  registra¬ 
tion  and  the  mapping  are  reliable  when  the  bandwidth  is 
larger  than  4  mm.  Computationally,  the  proposed  approach 
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Fig.  8.  Displacement  of  region  of  interest  boundary  points  as  a  function  of  respiration  phase  for  three  thoracic  patients,  (a) 
Mean  displacement  vs.  phase,  (b)  Maximum  displacement  vs.  phase.  4D  CT  =  four-dimensional  computed  tomography. 


resides  between  a  deformable  model-based  mapping  and 
a  surface  model-based  ROI  contour  mapping. 

The  success  of  the  image  content-based  approaches,  such 
as  the  proposed  narrow-band  approach  or  conventional  de¬ 
formable  image  registration,  arises  from  the  fact  that  they 
fully  utilize  the  inherent  image  features  of  the  two  input  im¬ 
ages.  The  narrow  band-based  technique  is  particularly  attrac¬ 
tive  because  it  takes  advantages  of  the  useful  features  of  both 
image  content-based  technique  and  the  regional  surface- 
based  model.  In  a  sense,  it  is  a  hybrid  approach  of  the  two  dis¬ 
tinct  types  of  algorithms.  The  narrow-band  approach  utilizes 
the  imaging  features  surrounding  the  ROI  to  guide  the  search 
of  the  optimal  mapped  contours  while  considering  the  shape 
integrity  of  the  ROI  surface.  It  eliminates  the  need  for  a  global 


registration  of  the  input  images  and  thus  greatly  increases  the 
computational  efficiency. 

Application  of  the  proposed  contour  mapping  technique  to 
five  clinical  cases  indicates  that  the  technique  is  accurate  and 
computationally  efficient.  A  common  problem  in  image 
segmentation  and  contour  mapping  studies  is  the  lack  of 
quantitative  validation.  In  the  studies  of  Lu  et  al.  (13)  and 
Schriebmann  et  al.  (14),  for  example,  the  accuracy  of 
a  deformable  model-based  contour  mapping  technique  was 
evaluated  purely  on  the  basis  of  visual  inspection.  Although 
it  is  a  convenient  way  for  rapid  assessment  of  a  segmentation 
calculation,  especially  in  a  case  in  which  the  “ground  truth” 
contours  do  not  exist,  the  method  falls  short  in  quantization. 
The  same  approach  was  used  in  many  other  previous 


Fig.  9.  Contour  propagation  in  a  head-and-neck  case,  (a)  Planning  computed  tomography  with  manually  outlined  template 
contours  (in  blue)  for  body,  mandible,  and  gross  tumor  volume,  (b)  Cone-beam  computed  tomography  along  with  contours 
after  warping  (in  red)  for  the  corresponding  structures. 
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investigations  (1,  5, 14,  30).  In  this  study,  a  bidirectional  con¬ 
tour  mapping  was  proposed  to  examine  the  reliability  and  ro¬ 
bustness  of  a  contour  mapping  technique.  This  method 
provides  a  useful  test  in  assessing  the  success  of  a  contour 
propagation  algorithm.  We  would  like  to  point  out  that  the  bi¬ 
directional  mapping  technique  introduced  in  this  work  is 
a  necessary  (but  not  sufficient)  test.  In  a  rare  but  possible  sit¬ 
uation,  the  bidirectional  mapping  may  not  be  able  to  find  that 
an  error  occurred  in  the  narrow-band  mapping  process.  A  vi¬ 
sual  inspection  of  the  mapped  result  may  help  in  this  situa¬ 
tion.  On  the  basis  of  the  bidirectional  mapping  experiments 
and  visual  inspection  for  the  patient  studies,  we  conclude 
that  the  proposed  approach  can  perform  very  well  even  in 
the  presence  of  significant  deformations. 

In  our  calculation,  we  observed  that  the  regular  grid  of 
B  Spline  control  points  could  be  mapped  to  a  region  outside 
the  narrow  band.  Although  it  seems  that  this  does  not  directly 
affect  the  accuracy  of  the  method,  it  may  prolong  the  calcu¬ 
lation  by  computing  the  displacements  in  regions  where  met¬ 
ric  information  is  irrelevant.  Setups  have  been  proposed  to 
adapt  the  splines  control  mesh  to  regions  where  deformation 
is  found  to  be  significant  (31),  and  the  extension  of  the 
method  would  allow  us  to  use  the  B  Spline  control  points  de¬ 
fined  only  in  the  regions  within  the  narrow  band.  Implemen¬ 
tation  of  this  type  of  technique  should  further  reduce  the 
computation  time  required  to  find  the  optimal  solution. 

Although  there  are  numerous  deformable  algorithms,  in¬ 
cluding,  for  example,  the  elastic  model  (32-34),  viscous  fluid 
model  (35),  optical  flow  model  (5,30,36),  finite  element 
model  (33,  37),  and  radial  basis  function  models  such  as 
the  basis  spline  model  (28,  38,  39)  and  thin  plate  spline  model 
(40-43),  a  truly  robust  tool  suitable  for  routine  clinical  appli¬ 
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cations  is  yet  to  be  developed.  Each  of  these  approaches  has 
its  pros  and  cons.  The  deformable  calculation  can  be  greatly 
facilitated  if  some  a  priori  system  information  can  be  incor¬ 
porated.  Along  this  line,  the  homologous  correspondence  of 
the  bony  structure  in  two  input  images  has  been  incorporated 
in  thin  plate  spline  method,  and  remarkable  improvement  has 
resulted  (44).  The  narrow  band-generated  ROI  contour  cor¬ 
respondence  could  also  be  used  as  prior  knowledge  to 
improve  a  deformable  registration.  This  work  is  still  in  prog¬ 
ress  and  will  be  reported  in  the  future. 

CONCLUSIONS 

In  this  work  we  have  developed  a  regional  deformable 
registration-based  method  to  auto-propagate  contours  for 
4D  radiotherapy.  The  central  idea  is  that  a  narrow  band 
encompassing  an  ROI  surface  carries  the  neighborhood  in¬ 
formation  of  the  ROI  surface  and  can  be  used  to  establish 
a  reliable  association  between  the  ROIs  in  two  phase-specific 
image  sets.  Different  from  other  type  of  regional  algorithms, 
such  as  surface  mapping,  the  method  uses  the  image  features 
captured  in  a  band  to  guide  the  search  for  the  optimal  contour 
mapping.  Compared  with  conventional  deformable  image 
registration-based  approaches,  a  great  reduction  in  computa¬ 
tional  burden  and  a  large  capture  radius  in  optimization  space 
result.  Our  study  demonstrated  that  the  information  contained 
in  the  boundary  region  can  be  used  to  guide  the  contour  map¬ 
ping  in  all  the  testing  cases  presented  in  this  article.  The  pro¬ 
posed  regional  model  decreases  the  workload  involved  in 
4D-CT  ROI  segmentation  and  provides  a  valuable  tool  for 
the  efficient  use  of  available  spatial-temporal  information 
for  4D  simulation  and  treatment  planning. 
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The  purpose  of  this  work  is  to  develop  a  novel  feature-based  registration  strategy  to  automatically 
map  the  rectal  contours  from  planning  computed  tomography  (CT)  (pCT)  to  cone  beam  CT 
(CBCT).  The  rectal  contours  were  manually  outlined  on  the  pCT.  A  narrow  band  with  the  outlined 
contour  as  its  interior  surface  was  then  constructed,  so  that  we  can  exclude  the  volume  inside  the 
rectum  in  the  registration  process.  The  corresponding  contour  in  the  CBCT  was  found  by  using  a 
feature-based  registration  algorithm,  which  consists  of  two  steps:  (1)  automatically  searching  for 
control  points  in  the  pCT  and  CBCT  based  on  the  features  of  the  surrounding  tissue  and  matching 
the  homologous  control  points  using  the  scale  invariance  feature  transformation;  and  (2)  using  the 
control  points  for  a  thin  plate  spline  transformation  to  warp  the  narrow  band  and  mapping  the 
corresponding  contours  from  pCT  to  CBCT.  The  proposed  contour  propagation  technique  is  applied 
to  digital  phantoms  and  clinical  cases  and,  in  all  cases,  the  contour  mapping  results  are  found  to  be 
clinically  acceptable.  For  clinical  cases,  the  method  yielded  satisfactory  results  even  when  there 
were  significant  rectal  content  changes  between  the  pCT  and  CBCT  scans.  Asa  consequence,  the 
accordance  between  the  rectal  volumes  after  deformable  registration  and  the  manually  segmented 
rectum  was  found  to  be  more  than  90%.  The  proposed  technique  provides  a  powerful  tool  for 
adaptive  radiotherapy  of  prostate,  rectal,  and  gynecological  cancers  in  the  future.  ©  2008  Ameri¬ 
can  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.2975230] 

Key  words:  image  guided  radiation  therapy  (IGRT),  image  registration,  deformable  model, 
segmentation,  scale  invariance  feature  transformation  (SIFT) 


I.  INTRODUCTION 

Patients  treated  with  radiotherapy  for  cancers  such  as  pros¬ 
tate,  rectal,  and  gynecological  cancers  experience  large  day- 
to-day  changes  in  their  rectal  volumes  due  to  motion,  disten¬ 
tion,  and  filling.  Due  to  variations  in  the  image  content,  an 
exact  correspondence  between  two  image  sets  acquired  at 
different  time  points  may  not  exist.  Thus,  any  deformable 
model  relying  on  the  use  of  information  contained  in  the 
entire  image  may  not  be  adequate  in  dealing  with  these  pa¬ 
tients.  The  artifacts- induced  disjoint  between  the  images  also 
makes  the  autopropagation  of  contours  outlined  in  one  set  of 
images  to  another  highly  difficult  with  conventional  strate¬ 
gies.  With  continued  enthusiasm  for  adaptive  radiotherapy, 
the  ability  to  reliably  and  efficiently  map  the  rectum  outlined 
in  the  planning  computed  tomography  (CT)  (pCT)  to  the 
on-treatment  cone  beam  CT  (CBCT)  images  now  becomes  a 
bottleneck  and  needs  to  be  resolved  in  order  for  many  pa¬ 
tients  with  cancer  within  the  pelvis  to  benefit  from  the  novel 
adaptive  replanning  strategy.1,2 

The  issue  of  rectal  motion  and  deformation  in  conformal 
radiation  therapy  is  described  in  various  publications.  Lee  et 
al.  evaluated  the  CBCT  as  a  tool  to  quantify  the  accuracy  and 
precision  of  a  simulated  IM  RT  treatment  delivery  model  for 
rectal  cancer  when  rectal  motion  due  to  filling  and  deforma¬ 
tion  was  taken  into  account.3  The  mean  deformation  varia¬ 
tion  of  0.71  and  0.94  cm  in  the  LAT  and  AP  directions  was 


reported.  Foskey  et  al.  shrank  the  rectal  gas  region  to  a  vir¬ 
tual  point  in  order  to  make  the  correspondence  of  the  rectal 
volumes  in  two  sets  of  images.4  Gao  etal.  used  an  automatic 
image  intensity  modification  procedure  to  create  artificial  gas 
pockets  in  the  pCT  images.5  The  major  drawbacks  of  these 
types  of  approaches  are  the  artificial  introduction  of  image 
features  within  the  rectal  volume  and  the  potentially  inaccu¬ 
rate  association  of  the  artificial  image  features.  As  a  conse¬ 
quence,  the  accordance  between  the  rectal  volumes  after  de¬ 
formable  registration  and  the  manually  segmented  rectum 
was  found  to  be  less  than  80%. 

In  this  work,  we  propose  to  use  the  image  information  in 
the  neighborhood  outside  the  rectal  wall  as  the  driving  force 
to  guide  the  rectal  contour  propagation  from  the  pCT  to 
CBCT.  Because  the  content  in  the  region  outside  the  rectal 
wall  should  be  conserved,  regardless  of  any  changes  in  the 
rectal  filling  and  distension,  this  strategy  seems  to  be  physi¬ 
cally  sensible.  Coupled  with  a  powerful  feature-based  de¬ 
formable  registration  model,  which  identifies  homologous 
tissue  features  shared  by  the  pCT  and  CBCT  images,  the 
novel  approach  captures  the  key  issues  of  the  system  and 
provides  a  natural  solution  to  the  above  stated  problem.  Ap¬ 
plication  of  the  proposed  algorithm  to  a  number  of  digital 
phantoms  and  clinical  cases  demonstrates  that  the  technique 
is  accurate  and  robust  and  may  be  useful  for  future  adaptive 
therapy  planning. 
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Fig.  1.  Overall  process  of  rectal  contour  propagation. 


II.  METHODS  AND  MATERIALS 
II.A.  Software  platform 

The  proposed  contour  mapping  algorithm  was  imple¬ 
mented  using  the  I nsight  Toolkit6  and  the  Visualization  Tool¬ 
kit  (VTK),7  which  are  open  source  cross-platform  C  +  +  soft¬ 
ware  toolkits  sponsored  by  the  National  Library  of  M  edicine. 
They  are  freely  available  for  research  purposes  (see  Refs.  34 
and  35).  ITK  provides  various  basic  algorithms  to  perform 
registration  and  segmentation  for  medical  images.  The  pro¬ 
grams  contained  in  ITK  are  highly  extendable,  making  it  an 
ideal  platform  for  development  of  image  registration  and 
processing  techniques.  VTK  is  primarily  used  for  image  vi¬ 
sualization  (including  contours). 

II.B.  Narrow  band  construction 

Inconsistency  in  rectal  contents  between  two  input  image 
sets  could  severely  reduce  the  performance  of  a  deformable 
registration  algorithm.  Coregistering  an  empty  rectum  with¬ 
out  bowel  gas  to  a  rectum  filled  with  bowel  gas  using  any 
deformable  model  could  be  problematic,  for  example.  A 
natural  strategy  is  to  exclude  the  volume  inside  the  rectal 
wall.  In  practice,  the  template  rectal  contour  in  the  pCT  im¬ 
age  has  been  manually  contoured  as  a  part  of  the  routine 
treatment  planning  process,  thus  making  it  a  straightforward 
matter  to  exclude  the  volume  inside  the  rectal  wall.  Figure  1 
shows  the  proposed  contour  mapping  process.  After  manual 
segmentation  on  the  pCT,  a  narrow  band  as  sketched  in  Fig. 
2  is  constructed  with  the  manually  segmented  rectum  repre¬ 
senting  the  inner  surface  of  the  band.  On  an  axial  slice,  the 
contour  has  a  polygon  shape  and  the  vertices  of  the  polygon 
form  the  basis  for  constructing  the  narrow  band.  The  distance 
between  the  neighboring  vertices  on  the  contour  is  typically 
2-10  mm  depending  on  the  shape  of  the  contour.  In  gener¬ 
ating  the  narrow  band,  we  first  create  squares  with  side 
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Fig.  2.  A  sketch  of  narrow  band,  (a)  A  narrow  band  image  surrounding  a 
manually  segmented  rectal  contour  and  (b)  a  narrow  band  construction  is 
illustrated  for  two  vertex  points  A  and  B. 


length  of  d  for  each  vertex,  as  depicted  by  points  A  and  B  in 
Fig.  2(b).  In  order  to  obtain  asmooth  band,  betweenA  and  B 
three  more  squares,  cornered  at  points  C,  D,  and  E,  are  in¬ 
serted.  Point  C  is  chosen  to  be  the  middle  point  between  A 
and  B.  Point  D  is  the  point  between  A  and  C,  and  point  E  is 
the  point  between  B  and  C.  More  interpolated  vertex  points 
can  be  similarly  introduced  to  obtain  a  smooth  band.  The 
principle  of  the  narrow  band  diameter  selection  is  to  exclude 
most  bony  structures  outside  the  narrow  band,  since  the  bony 
structures  are  rigid  and  heavily  affect  the  control  point  selec¬ 
tion.  M  eanwhile,  the  generated  narrow  band  can  capture  suf¬ 
ficient  information  to  drive  the  finding  of  its  counterpart  in 
the  subsequent  CBCT.  In  general,  the  size  of  the  squares  is 
therefore  within  1  cm,  so  that  the  diameter  of  the  narrow 
band  is  within  1.5  cm. 

The  narrow  band  in  our  approach  is  used  as  a  compact 
representation  of  the  rectal  surface.  As  will  be  detailed  in  the 
next  subsection,  a  feature-based  deformable  registration  al¬ 
gorithm  is  employed  to  find  the  correspondence  of  the  band 
in  theCBCT  images.  Upon  successful  registration,  the  defor¬ 
mation  field  is  utilized  to  propagate  the  pCT  contour  to  the 
CBCT.  Because  only  the  image  features  outside  the  rectum 
are  used,  a  narrow  band  shown  in  Fig.  2  permits  us  to  take 
advantage  of  the  regional  information  inside  the  narrow  band 
yet  avoiding  the  nuisance  of  rectum/bladder  filling. 

II.C.  Feature-based  warping  of  the  narrow  band 

As  illustrated  in  Fig.  1,  the  process  of  contour  mapping  is 
to  warp  the  narrow  band  constructed  above  in  such  a  way 
that  its  best  match  in  the  CBCT  images  is  found.  M  athemati- 
cally,  this  constitutes  an  optimization  problem,  in  which  a 
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group  of  transformation  parameters  transform  the  points 
within  the  band  in  the  pCT  to  their  corresponding  points  in 
the  CBCT.  The  input  to  the  contour  mapping  software  in¬ 
cludes  the  narrow  band  and  the  CBCT  images,  which  are 
described  by  the  image  intensity  distributions  la(x)  and 
lb(x),  respectively. 

To  find  the  transformation  matrix,  T(x),  that  maps  an  ar¬ 
bitrary  point  in  the  band  to  the  corresponding  point  in  the 
CBCT  images  (or  vice  versa),  a  thin  plate  spline  (TPS)  de¬ 
formable  model  is  employed.  But  other  models  should  also 
be  applicable  to  model  the  deformation  of  the  band.  We  au¬ 
tomates  the  control  point  selection  by  using  the  scale  invari¬ 
ance  feature  transformation  (SIFT)  tissue  feature  searching 
(see  next  subsection  for  details).  Roughly,  300  control  points 
are  selected  based  on  the  prominent  tissue  features. 

The  detailed  description  of  the  TPS  transformation  can  be 
found  in  Ref.  8.  For  two-dimensional  (2D)  images,  a  weight¬ 
ing  vector  W  =(w1,w2, ...  ,w„)  and  the  coefficients  ai,au,av 
are  computed  from  a  series  of  matrices  which  are  constructed 
using  n  pairs  of  selected  control  points  in  the  fixed  image 
(x/(y()  and  in  the  moving  image  (u, , f,),  respectively.  The 
function  transforming  a  pixel  coordinate  in  the  moving  im¬ 
age  to  a  new  coordinate  in  the  fixed  image  is  defined  as 

n 

f(u',v')  =  a:  -Fa^  +avv  +2  (|p,  - (u,t-)|),  (1) 

i=  o 

where  pf  is  the  control  points  coordinate  in  the  fixed  image 
and  U  is  a  basis  function  to  measure  the  distance. 


II.D.  SIFT 

The  feature-based  deformable  registration  is  an  essential 
part  of  the  proposed  contour  mapping  process.  FI  ere,  we  au¬ 
tomate  the  control  point  selection  by  using  the  SIFT-based 
tissue  feature  searching.  Because  of  the  efficient  use  of 
a  priori  system  knowledge,  the  approach  greatly  enhances 
the  robustness  of  the  narrow  band  warping  algorithm. 

The  SIFT  method  was  introduced  by  Lowe  to  characterize 
the  local  tissue  features.  The  method  utilizes  both  image  in¬ 
tensity  and  local  gradient  information  to  characterize  the 
neighborhood  property  of  a  point.9  The  algorithm  includes 
scale-space  extrema  detection,  control  point  localization,  ori¬ 
entation  assignment,  and  control  point  descriptor.  In  2D 
cases,  for  example,  the  method  uses  the  orientation  histo¬ 
grams  of  the  four  quadrants  surrounding  a  point  (containing 
64  pixels)  to  characterize  the  inherent  tissue  feature  of  the 
point  (see  Fig.  3).  To  obtain  the  histogram  for  a  quadrant,  as 
illustrated  in  Fig.  3,  the  gradient  of  each  of  the  16  pixels  in  a 
quadrant  is  computed.  An  eight-bin  histogram,  with  first  bin 
representing  the  number  of  pixels  whose  gradients  fall  be¬ 
tween  0°  and  45°,  and  so  forth,  is  then  constructed.  For  il¬ 
lustration,  the  histogram  of  each  of  the  four  quadrants  is 
displayed  schematically  in  the  right  panel  of  Fig.  3  as  an 
eight-vector  plot.  In  total,  32  vectors  are  calculated  in  2D 
case.  In  extending  the  SIFT  method  from  2D  to  three  dimen¬ 
sional  (3D),  total  of  192  vectors  are  needed.  These  vectors 
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Fig.  3.  A  sketch  of  orientation  histogram  in  SIFT  method.  The  gradient  of 
each  of  the  16  pixels  in  a  quadrant  is  computed.  An  eight-bin  histogram, 
with  first  bin  representing  the  number  of  pixels  whose  gradients  fall  between 
0°  and  45°,  and  so  forth,  is  then  constructed.  The  histogram  of  each  of  the 
four  quadrants  is  displayed  schematically  in  the  right  panel  as  an  eight- 
vector  plot. 


represent  the  local  feature  and  serve  as  a  signature  of  the 
point.  The  SIFT  descriptor  is  considered  as  one  of  the  most 
effective  descriptors  currently  available.10,11 

Theoretically,  the  SIFT  descriptor  can  be  computed  for 
each  voxel  in  an  image.  However,  this  is  computationally 
expensive.  The  commonly  used  sampling  strategy  is  to  com¬ 
pute  the  descriptor  every  2-3  voxels  in  x,  y,  and  z  directions. 
After  the  SIFT  descriptors  are  computed  in  both  input  im¬ 
ages,  the  points  having  the  most  similar  SIFT  descriptors  in 
the  two  images  are  then  identified.  For  a  given  point,  indexed 
by  n,  in  the  pCT  image,  the  least-squares  difference  of  the 
SIFT  descriptor  of  the  point  and  that  of  a  potential  associa¬ 
tion  point  n’  in  the  CBCT,  SM.,  is  first  computed  according 
to 

Sn,n'  =  A/E|(V/n)a-(V/n,)J2,  (2) 

where  /  represents  the  image  intensity,  a  indexes  the  bins  of 
the  SIFT  histogram  of  a  point  and  the  summation  over  a  runs 
from  1  to  32  for  the  2D  case,  and  1  to  192  for  the  3D  case. 
Typically,  about  1000  SIFT  descriptors  n,n’  are  computed  in 
the  narrow  band  in  the  pCT  and  CBCT,  respectively.  It  is 
unnecessary  to  determine  SM,  for  all  possible  combinations 
n,n',  which  may  dramatically  increase  the  calculation  time. 
We  use  a  specific  search  radius  to  control  the  number  of  S„x 
calculation.  The  mapping  results  are  more  accurate  with 
larger  search  radius,  however,  the  calculation  time  of  SIFT 
mapping  becomes  longer.  After  Snn.  is  computed,  two  points 
n[  and  n2  that  have  the  least  histogram  difference  with  point 
n  are  identified.  If  the  ratio  (for  convenience,  the  ratio  is 
referred  to  as  the  k  ratio  hereafter)  of  these  two  values  is  less 
than  80%,  the  point  that  has  the  least  5  value  is  chosen  ten¬ 
tatively  as  the  correspondence  of  the  point  n,  otherwise,  no 
association  is  made  for  the  point.  The  k  ratio  varies  between 
0  and  1  and  is  an  empirical  measure  of  feature  correspon¬ 
dence  between  two  images.  The  lower  the  k  ratio,  the  "stron¬ 
ger"  the  association  of  the  two  feature  points  on  pCT  and 
CBCT.  Because  of  the  inherent  difference  in  the  textures  of 
the  involved  organs,  the  determination  of  the  k  ratio  may  be 
organ  specific.  Typically,  it  is  determined  by  a  tradeoff  be- 
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tween  the  number  of  associated  point  pairs  and  the  reliability 
of  the  associations.  For  a  bone,  the  feature  is  clear  and  hun¬ 
dreds  point  pairs  can  be  associated  under  a  threshold  of  50%. 
On  the  other  hand,  for  the  rectum,  the  feature  is  not  as  ob¬ 
vious  as  bone.  If  we  still  use  this  low  k  ratio,  the  number  of 
association  pairs  may  be  very  limited.  In  this  situation,  a 
higher  threshold,  say  80%,  is  usually  used  to  increase  the 
number  of  associated  point  pairs. 

To  further  increase  the  accuracy  of  feature  point  associa¬ 
tion,  a  bidirectional  mapping  strategy  is  developed  based  on 
the  fact  that  if  a  point  in  the  pCT  is  mapped  correctly  to  the 
CBCT,  it  will  be  default  to  be  mapped  back  to  the  original 
point  in  the  pCT  when  an  inverse  map  is  applied  to  the 
corresponding  point  in  the  CBCT.  Therefore,  after  the  origi¬ 
nal  association  of  feature  points  as  described  above,  the 
mapped  points  in  CBCT  are  inversely  coregistered  to  the 
pCT.  If  the  correspondence  still  exists,  the  associated  point 
pair  is  labeled  a  match.  Otherwise,  they  are  considered  as  a 
mismatch  and  deleted  from  the  list  of  correspondence  points. 
Upon  the  association  of  the  feature  points,  the  associated 
points  are  employed  as  control  points.  The  control  point  in 
pCT  and  CBCT  corresponds  each  other,  thus  the  numbers  of 
control  points  in  the  two  input  images  are  the  same.  It  was 
noticed  that,  when  the  CBCT  region  of  interest  (ROI)  is  ex¬ 
panded,  the  increase  of  feature  point  generation  does  not 
affect  the  control  point  association  and  final  contour  map¬ 
ping.  The  coordinates  of  an  arbitrary  point  on  the  contour  in 
CBCT  are  obtained  by  interpolating  the  displacement  vectors 
of  the  control  points  using  TPS  transformation  after  the  con¬ 
trol  point  association  is  established. 


II.E.  Evaluation  of  the  models  using  digital  phantom 
and  existing  patient  data 

The  performance  of  the  above  model  is  evaluated  by  a 
number  of  2D  digital  phantoms  and  archived  clinical  cases. 
In  the  digital  phantom  experiments,  two  deformations  are 
introduced.  A  virtue  of  this  approach  is  that  the  "ground 
truth"  solutions  exist  and  the  transformation  matrices  are 
known,  thus  making  the  evaluation  straightforward.  The 
mathematical  transformations  used  to  deform  the  phantom 
are  generated  using  a  formula12 

x'(x,y)=(l+ibcosm0)x,  (3) 


y'(x,y)  =(1  +b  cos  m0)y.  (4) 

Here,  6>=tan-1  y/x.  Two  parameters,  m  and  b,  are  used  to 
characterize  a  deformation.  Generally,  they  describe  the 
complexity  and  magnitude  of  a  deformation,  respectively. 
The  contour  outlined  in  the  original  image  is  then  mapped  to 
the  deformed  image.  The  accuracy  of  the  contour  mapping 
calculation  is  assessed  by  comparing  directly  with  the  de¬ 
formable  mapping  from  the  known  transformation  matrix. 

Contour  propagation  from  pCT  to  CBCT  is  studied  by 
using  three  prostate  cancer  patients  and  two  rectal  cancer 
cases.  The  pCT  is  acquired  with  a  GE  Discovery-ST  CT 
scanner  (GE  Medical  System,  Milwaukee,  Wl)  approxi¬ 
mately  two  weeks  prior  to  the  initiation  of  the  radiotherapy. 


The  on-treatment  CBCT  images  are  acquired  using  the 
Varian  Trilogy™  (Varian  Medical  Systems,  Palo  Alto,  CA). 
Each  slice  of  pCT  or  CBCT  is  discretized  into  512x512 
voxels.  The  images  are  transferred  through  DICOM  to  a 
high-performance  personal  computer  with  a  Xeon  (3.6  GHz) 
processor  for  image  processing.  The  manually  outlined  con¬ 
tours  in  the  pCT  images  are  mapped  to  CBCT  images  using 
the  proposed  technique.  For  the  cases  studied  here,  the  ac¬ 
cordance  between  the  rectal  volumes  after  deformable  regis¬ 
tration  and  the  manually  segmented  rectum  is  employed  to 
assess  the  success  of  the  proposed  algorithm. 

To  quantitatively  evaluate  the  result  of  contour  propaga¬ 
tion,  the  accordance  value  between  the  automapped  contour 
and  manually  outlined  contours  were  calculated.  In  general, 
suppose  A  and  B  are  two  contours,  the  accordance  value  r  is 
defined  as 


Va  n  V_b 

vauvb' 

here,  V  is  the  containing  volume  of  A  or  B. 


(5) 


III.  RESULTS 

III.A.  2D  digital  phantom  experiment 

The  proposed  algorithm  is  first  tested  using  a  2D  digital 
phantom  [Fig.  4(a)]  with  two  intentionally  introduced  defor¬ 
mations  of  the  image  shown  in  Figs.  4(b)  and  4(c),  respec¬ 
tively.  The  rectal  contour  is  manually  outlined  and  shown  in 
Fig.  4(a).  The  deformation  shown  in  Figs.  4(b)  and  4(c)  are 
obtained  by  setting  the  parameters  b  and  m  in  Eqs.  (3)  and 
(4)  to  (b=  2,  m=  2)  and  (b=  2,  m=  3),  respectively.  The  curves 
close  to  the  interior  surface  of  the  rectum  in  Figs.  4(b)  and 
4(c)  represent  the  automapped  contour.  For  comparison,  the 
original  contour  in  Fig.  4(a)  is  also  mapped  rigidly  to  Figs. 
4(b)  and  4(c).  Overall,  the  mapped  contours  can  capture  the 
main  features  of  the  two  dramatic  deformations,  and  conform 
to  the  boundary  of  the  rectum  in  both  cases. 

In  obtaining  the  result  shown  in  Fig.  4(b),  a  total  of  200 
control  points  were  identified  by  the  bidirectional  SIFT  cal¬ 
culation  as  described  in  method.  Note  that  the  bony  structure 
in  the  image  has  been  excluded  in  this  calculation  by  setting 
an  intensity  threshold  of  300  CT  number.  In  this  way,  any 
unphysical  bony  structure  deformation  is  avoided.  For  clar¬ 
ity,  a  selection  of  the  SI FT-identified  control  point  associa¬ 
tions  are  displayed  in  Fig.  5.  The  superior  contour  represents 
the  superior  surface  of  narrow  band.  The  total  number  of 
control  points  identified  here  are  far  more  than  that  com¬ 
monly  used  in  TPS  calculation,13  allowing  an  improved  de¬ 
formable  warping  of  the  narrow  band.  We  should  notice  that 
the  control  points  2,  3,  4,  5,  and  9  in  Fig.  5  are  relatively  far 
away  from  the  rectum  wall  compared  to  control  points  6  and 
10.  Since  the  TPS  interpolation  is  used  after  SIFT  mapping, 
every  control  point  including  points  2,  3,  4,  5,  and  9  will 
affect  the  deformable  warping  and  therefore  the  contour 
shape,  although  the  weights  of  points  2,  3,  4,  5,  and  9  are 
smaller  than  points  6  and  10.  The  displacement  field  derived 
by  using  TPS  method  is  shown  in  Fig.  6(a).  For  comparison, 
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(c) 


Fig.  4.  Rectal  contour  propagation  from  the  2D  pCT  slice  to  two  dramati¬ 
cally  deformed  images,  (a)  Original  contour,  (b)  and;  (c)  its  optimal  map¬ 
ping  in  the  two  deformed  images. 


the  known  displacement  field  from  Eqs.  (3)  and  (4)  is  plotted 
in  Fig.  6(b).  The  subtraction  between  the  TPS-derived  dis¬ 
placement  field  and  the  known  field  is  shown  in  Fig.  6(c).  It 
is  found  that  the  average  deviation  of  the  SIFT-TPS  displace¬ 
ment  from  the  known  solution  is  less  than  1.2  mm. 

III.B.  Clinical  case  study 

The  contour  propagation  study  from  pCT  to  CBCT  for  the 
first  prostate  case  is  presented  in  Fig.  7.  The  top  row  shows 
the  pCT  image  with  manually  outlined  contours.  The  au- 
tomapped  contours  overlaid  on  the  CBCT  are  displayed  in 
the  bottom  row.  For  comparison,  the  manually  outlined  con¬ 
tours  on  the  CBCT  are  also  plotted  in  the  bottom  row.  As 
mentioned  in  Sec.  I,  the  propagation  of  rectum  wall  is  often 
complicated  by  the  fact  that  the  physical  one-to-one  corre¬ 
spondence  may  not  exist  due  to  the  addition  or  subtraction  of 
some  contents  within  the  rectum.  Figure  8  exemplifies  this 
and  shows  that  the  rectal  filling  at  the  time  of  CBCT  acqui¬ 
sition  is  quite  different  from  that  of  pCT.  As  can  be  intu¬ 
itively  conceived,  this  image  content  change  could  severely 


(b) 


Fig.  5.  Control  points  in  the  2D  contour  mapping. 


reduce  the  performance  of  a  conventional  deformable 
registration.14  The  narrow  band  approach  described  in  this 
work  circumvents  the  problem  by  excluding  the  rectal  vol¬ 
ume  affected  by  the  rectum/bladder  filling.  Accuracy  was 
evaluated  by  comparison  with  manually  outlined  contours  on 
the  CBCTs. 15-17  It  is  clear  that  the  mapped  contours  closely 
conform  to  the  rectal  wall  change.  The  accordance  between 
the  rectal  volume  extended  by  the  automapped  contour  and 
the  manually  segmented  rectal  volume  was  found  to  be  more 
than  90%. 

In  practice,  rectal  volume  motion  and  deformation  can 
cause  large  uncertainties  pertaining  to  the  adequacy  of  actual 
dose  delivered  to  the  gross  tumor  volume  as  well  as  to  the 
surrounding  normal  structures.  This  issue  has  been  a  major 
obstacle  in  the  implementation  of  IMRT  in  rectal  cancer.  In 
Fig.  8  six  axial  pCT  and  CBCT  images  of  a  rectal  cancer 
patient  acquired  in  an  interval  of  two  weeks  are  shown. 
Large  target  volume  motion  and  deformation  are  observed 
from  Fig.  8.  The  rectal  volume  in  the  pCT  is  found  to  be 
three  times  more  than  that  of  the  rectal  volume  in  the  CBCT 
and  thus  represents  a  challenging  situation  for  any  deform¬ 
able  model.  The  rectal  contours  are  manually  outlined  in  the 
pCT  and  mapped  to  the  subsequent  CBCT  using  the  pro¬ 
posed  method.  The  first  and  second  rows  of  Fig.  8  show  six 
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Fig.  6.  Displacement  fields,  (a)  TPS-derived  displacement  field  for  the  2D 
digital  phantom  study;  (b)  intentionally  introduced  displacement  field;  and 
(c)  subtraction  of  TPS  derived  and  the  known  displacement  fields. 


axial  slices  of  the  pCT  with  manually  outlined  contours.  The 
results  of  contour  propagation  from  the  pCT  to  the  CBCT  are 
shown  in  the  third  and  fourth  rows  of  Fig.  8.  As  the  same  as 
in  Fig.  7  the  manually  outlined  contours  on  the  CBCT  are 
also  plotted.  The  accordance  between  the  rectal  volume  ex¬ 
tended  by  the  automapped  contours  and  the  manually  seg¬ 
mented  rectal  volume  was  found  to  be  more  than  95%.  The 
rectal  deformations  in  Fig.  8  are  quite  large  and  thus  present 
challenges  to  any  deformable  model  or  contour  mapping 
technique.  It  is  impressive  that  a  simple  approach  with  a 
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narrow  band  and  SIFT  descriptor  can  capture  the  main  fea¬ 
ture  of  the  rectal  contour  and  help  to  find  the  correspondence 
contours  in  the  CBCT  images. 

To  further  examine  the  performance  of  the  proposed  tech¬ 
nique,  the  method  was  also  applied  to  three  additional  pa¬ 
tients  (Fig.  9).  The  automapped  contours  are  plotted  together 
with  the  manually  outlined  contour  on  the  CBCT.  For  com¬ 
parison,  the  original  contours  on  pCT  are  also  mapped  rig¬ 
idly  to  the  CBCT.  The  accordance  values  between  the  pCT 
and  CBCT  contours,  as  well  as  between  the  automapped  and 
manually  segmented  CBCT  contours  for  these  three  patients 
are  listed  in  Table  I.  In  these  cases,  the  accordance  values  are 
increased  from  around  75%  to  over  90%  after  contour  map¬ 
ping.  The  influence  of  the  k  ratio  on  the  contour  propagation 
is  illustrated  by  the  data  listed  in  Table  II,  where  the  accor¬ 
dance  values  for  a  few  different  k  ratios  for  the  three  patients 
are  shown.  The  accordance  reached  its  peak  value  when  the 
k  ratio  is  between  0.8  and  0.9  for  all  these  three  cases.  When 
the  k  ratio  is  lower  than  0.8,  the  accordance  decreases  with 
the  decrease  of  the  k  ratio  because  less  control  points  are 
selected.  The  accordance  also  decreases  with  the  increase  of 
the  k  ratio  for  the  k  ratio  higher  than  0.9.  The  accordance 
value  is  stable  for  k  ratios  between  0.8  and  0.9.  The  data  also 
indicate  that  the  k  ratio  is  generally  organ  specific  and  is 
insensitive  for  different  patients. 


IV.  DISCUSSION 

In  this  work,  an  effective  feature-based  rectal  contour 
mapping  algorithm  has  been  described.  An  indispensable 
step  toward  online  or  offline  adaptive  replanning  with  con¬ 
sideration  of  the  patient's  dose  delivery  history  and  on- 
treatment  anatomy  is  the  expedite  organ  segmentation  of 
CBCT  images.  While  this  task  is,  in  principle,  achievable 
using  deformable  registration  of  the  pCT  and  CBCT  images, 
the  accuracy  of  the  registration  and  therefore  the  contour 
mapping,  is  often  adversely  affected  by  the  presence  of  im¬ 
age  contents  in  one  image  that  do  not  have  correspondence 
in  the  other  image.  The  propagation  of  rectum  wall  is  an 
example  of  this.  For  prostate,  rectal,  or  gynecological  cancer 
patients  for  example,  the  presence  and  absence  of  bowel  gas 
can  vary  daily.  Coregistering  an  empty  rectum  without  bowel 
gas  to  a  rectum  filled  with  bowel  gas  (or  vice  versa)  using 
any  deformable  model  could  be  problematic  and  large  errors 
could  occur. 

We  describe  a  regional  contour  propagation  algorithm  tak¬ 
ing  into  account  possible  organ  deformation  and  anatomic 
changes.  Because  the  narrow  band  contains  only  the  image 
features  outside  the  rectum,  this  method  is  not  affected  by 
the  rectum  filling  changes.  The  use  of  SIFT  descriptor  en¬ 
hances  our  ability  to  find  the  correspondence  of  the  narrow 
band  because  of  the  effective  utilization  of  image  intensity 
and  gradient  information.  In  contrast  to  the  conventional 
intensity- based  image  registration,  which  only  uses  intensity 
information  of  the  voxels,  the  feature-based  registration  ex¬ 
tracts  information  regarding  image  structure,  including 
shape,  texture,  etc.  Therefore,  the  feature-based  image  regis- 
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Fig.  7.  3D  contour  mapping  for  the 
rectum  of  a  man  with  prostate  cancer. 
The  top  row  is  the  three  transactions  in 
the  planning  CT  image,  the  bottom 
row  is  corresponding  transactions  in 
the  CBCT  image.  The  left  column  is 
the  axial  plane,  the  middle  column  is 
the  coronal  plane,  and  the  right  col¬ 
umn  is  the  sagittal  plane. 


tration  is  generally  more  effective  in  correctly  identifying 
corresponding  voxels  compared  to  the  intensity-based  image 
registration. 

In  this  study,  a  bidirectional  SIFT  descriptor  is  employed 
to  examine  the  reliability  and  robustness  of  the  calculations. 
The  bidirectional  mapping  further  enhances  the  degree  of 
success  of  a  contour  propagation  algorithm.  It  is  useful  to 
note  that  the  bidirectional  mapping  is  a  necessary  (but  not 


sufficient)  test.  In  a  rare  but  possible  situation,  the  bidirec¬ 
tional  mapping  may  not  be  able  to  find  an  error  occurred  in 
the  contour  mapping  process. 

Because  the  iterative  procedure  in  the  B-spline  is  not 
needed  in  our  method,  the  calculation  speed  is  at  least  ten 
times  faster  than  B-spline  registration.  Typically  the  total  cal¬ 
culation  time  of  SIFT-TPS  mapping  with  about  1000  SIFT 
descriptors  is  less  than  2  min.  Several  parameters  influence 


Fig.  8.  Rectal  contour  mapping  for  a 
rectal  cancer  case.  The  first  and  second 
rows  show  six  axial  slices  in  the  pCT 
image.  The  third  and  fourth  rows  are 
the  corresponding  slices  in  the  CBCT 
image. 
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(a) 


(b) 


(c) 


Fig.  9.  Automatic  contour  propaga¬ 
tion  for  three  additional  patients. 


(d) 


(e) 


(f) 


(g) 


(h) 


(i) 


calculation  time.  For  example,  larger  narrow-band  will  result 
in  longer  calculation  times.  The  number  of  control  points 
also  affects  calculation  time  quite  nonlinearly  as  well.  For 
most  cases,  300  control  points  are  enough  for  accurate  con¬ 
tour  mapping.  Due  to  the  tight  clinical  timeframes  (espe¬ 
cially  for  real-time  adaptive  schemes),  1  or  2  min  calculation 
time  allows  the  use  of  the  contour  mapping  tool  between 
acquiring  the  verification  images  and  delivering  the  dose 
fraction  for  online  corrections.18 

In  some  cases,  no  corresponding  feature  is  found  by  SIFT 
in  a  certain  area  close  to  the  rectal  wall.  For  instance,  no 
control  point  was  found  in  the  upper  part  of  Fig.  5.  Since  no 
large  local  deformation  was  in  these  regions,  the  result  was 
all  right.  However,  it  would  have  lead  to  larger  errors  in  case 


of  large  deformations  plus  low  feature  density,  which  may 
happen  in  smooth  soft  tissues.  We  will  improve  it  in  our 
future  work. 

One  of  the  practical  concerns  is  that  the  relatively  low 
quality  of  CBCT  images  may  influence  the  accuracy  of  im¬ 
age  registration  and  thus  the  contour  mapping.  Paquin  et  al. 
quantitatively  studied  the  influence  of  different  types  of 
noises  on  deformable  registration  and  found  that  the  accu¬ 
racy  of  image  registration  does  not  depend  on  the  global 
noise  unless  the  noise  reaches  a  certain  threshold  value.19 
Murphy  et  al.  also  demonstrated  that  noise  levels  in  cone- 
beam  CTs  that  might  reduce  manual  contouring  accuracy  do 
not  reduce  image  registration  and  automatic  contouring 
accuracy.20 


Table  I.  Accordance  values  between  the  pCT  and  CBCT  contours,  as  well  as  between  the  auto-mapped  and 
manually  segmented  CBCT  contours  for  three  patients. 


Patient  1 

Patient  2 

Patient  3 

Accordance  values  between  pCT  and 

CBCT  contours  (%) 

73.7 

76.5 

76.3 

Accordance  values  between  auto-mapped 
and  manually  segmented  CBCT  contours  (%) 

93.3 

91.3 

91.4 
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Table  II.  Accordance  values  between  the  auto-mapped  and  manually  segmented  CBCT  contours  for  different 
k-  rati  os  for  three  patients. 


k  ratio 

0.7 

0.75 

0.8 

0.85 

0.9 

0.95 

Patient  1  (%) 

81.5 

88.4 

91.6 

92.6 

93.3 

82.9 

Patient  2  (%) 

76.3 

82.1 

91.3 

88.9 

85.8 

79.1 

Patient  3  (%) 

67.5 

86.6 

88.3 

88.4 

91.4 

62.6 

Deformable  model  plays  an  important  role  in  automated 
contour  propagation.  N  umerous  approaches  have  been  devel¬ 
oped  for  different  applications.  Most  popular  deformable 
registration  methods  for  medical  images  include  the  thin 
plate  splines  (TPS),8  B -splines,21,22  and  finite  element 
method  (FEM  ).23  TPS  is  less  sensitive  to  noise  because  of  its 
global  calculation  nature.24  It  relies  on  the  use  of  homolo¬ 
gous  control  points  in  the  two  input  image  sets  to  be  coreg¬ 
istered.  Control  points  are  manually  selected  for  many  TPS 
applications.13,25,26  This  may  introduce  interuser  variability 
and  is  a  major  source  of  error.  Malsch  et  al.  presented  an 
automatic  block  matching  method,18  which  is  similar  to  the 
control  volume  based  approach  proposed  by  Schreibmann 
and  Xing.27  Kim  et  al.  presented  an  automated  TPS,  where 
an  arbitrary  set  of  control  points  is  supplied  initially  and  then 
is  iteratively  repositioned  until  the  resulting  warp  optimizes 
some  measure  of  registration.28"30  The  convergence  of  the 
iterative  calculation  is  slow  because  each  control  point  influ¬ 
ences  the  transformation  in  a  global  fashion.  An  alternative  is 
to  use  B-splines.  In  contrast  to  TPS,  which  allows  arbitrary 
configurations  of  the  control  points,  B -spline  requires  a  regu¬ 
lar  mesh  of  control  points  with  uniform  spacing.  Unlike 
spline-based  registration  methods,  FEM  models  the  deform¬ 
ing  image  as  an  elastic  body  subject  to  external  forces  which 
drive  the  deformation  and  internal  forces  (stresses)  which 
impose  smoothness  constraints.31,32  FEM  may  fail  to  model 
highly  localized  deformations,  since  the  deformation  energy 
caused  by  stress  increases  proportionally  with  the  strength  of 
the  deformation.33 


V.  CONCLUSION 

Large  interactional  patient  setup  uncertainty  and 
anatomy  changes  have  been  reported  in  numerous  studies 
and  are  widely  recognized  as  one  of  the  major  limiting  fac¬ 
tors  for  maximum  exploitation  of  modern  radiation  therapy 
techniques  such  as  IM  RT  and  IGRT.  The  advent  of  onboard 
volumetric  imaging  devices  promises  to  improve  the  situa¬ 
tion  by  providing  valuable  3D  (or  even  possibly  four¬ 
dimensional)  geometric  data  of  the  patient  in  the  treatment 
position  and  allows  for  the  adaptive  modification  of  treat¬ 
ment  plan  during  a  course  of  treatment. 

In  this  work,  an  effective  feature-based  rectal  contour 
mapping  algorithm  has  been  described.  The  method  yielded 
satisfactory  mapping  for  both  digital  phantom  and  clinical 
cases.  It  is  impressive  that  the  algorithm  is  able  to  success¬ 
fully  map  the  contours  from  pCT  to  CBCT  even  for  some 
very  challenging  cases  in  which  the  deformation  and/or  im¬ 


age  content  change  are  dramatic.  The  two  salient  features  of 
the  described  algorithm  are:  (1)  the  use  of  inherent  tissue 
feature  for  control  point  selection  as  a  priori  knowledge  for 
deformable  registration;  and  (2)  limiting  the  ROI  to  exclude 
the  volume  inside  the  rectum  and  focusing  on  the  adjacent 
neighborhood  of  the  rectal  contour.  The  algorithm  should  be 
extendable  for  contour  propagation  of  organs  with  similar 
features,  such  as  the  bladder  and  stomach. 

ACKNOWLEDGMENTS 

This  work  was  supported  in  part  by  grants  from  the 
Department  of  Defense  (W81XWH -06-1-0235  and 
W81XWH 05-1-0041)  and  National  Cancer  Institute  (1R01 
CA  98523  and  C A 104205). 

a)E lectronic  mail:  lei@reyes.stanford.edu 

1L.  Xing  eta/.,  "Overview  of  image-guided  radiation  therapy,"  Med.  Do- 
sim.  31  91-112  (2006). 

2A.  de  la  Zerda,  B.  Armbruster,  and  L.  Xing,  "Formulating  adaptive  ra¬ 
diation  therapy  (ART)  treatment  planning  into  a  closed-loop  control 
framework,"  Phys.  Med.  Biol.  5 2,  4137-4153  (2007). 

3P.  Lee etal.,  "Image-guided  radiation  therapy  (RT)  for  rectal  cancer  using 
cone  beam  CT  (CBCT),”  Int.  J.  Radiat.  Oncol.,  Biol.,  Phys.  66  S276 
(2006). 

4M .  Foskey  et  al.,  "Large  deformation  three-dimensional  image  registra¬ 
tion  in  image-guided  radiation  therapy,"  Phys.  Med.  Biol.  5Q  5869-5892 
(2005). 

5S.  Gao  et  al.,  "A  deformable  image  registration  method  to  handle  dis¬ 
tended  rectums  in  prostate  cancer  radiotherapy,"  Med.  Phys.  33  3304- 
3312  (2006). 

6L.  Ibanez,  W.  Schroeder,  and  L.  Ng,  "ITK  Software  Guide,”  Kitware, 
Inc.,  2003. 

7W.  Schroeder,  K.  M  artin,  and  B.  Lorensen,  The  Visualization  Toolkit:  An 
Object-Oriented  Approach  To  3D  Graphics,  3rd  ed.,  2003. 

8F.  L.  Bookstein,  "Principal  warps:  Thin  plate  splines  and  the  decomposi¬ 
tion  of  deformations,"  IEEE  Trans.  Pattern  Anal.  Mach.  Intell.  U  567- 
585  (1989). 

9D.  G.  Lowe,  "Object  recognition  from  local  scale-invariant  features," 
Proc.  of  the  International  Conference  on  Computer  Vision,  Corfu,  1999. 
10K.  Mikolajczyk  and  C.  Schmid,  "A  performance  evaluation  of  local  de¬ 
scriptors,"  IEEE  Trans.  Pattern  Anal.  Mach.  Intell.  27, 1615-1630  (2005). 
nG.  Wu,  F.  Qi,  and  D.  Shen,  "A  general  learning  framework  for  non-rigid 
image  registration,"  M  IAR,  2006,  pp.  219-227. 

12W.  Lu  etal.,  "Fastfree-form  deformable  registration  via  calculus  of  varia¬ 
tions,"  Phys.  Med.  Biol.  49,  3067-3087  (2004). 

13J.  Lian  et  al.,  "Mapping  of  the  prostate  in  endorectal  coil-based  MRI/ 
MRSI  and  CT:  A  deformable  registration  and  validation  study,"  Med. 
Phys.  31  3087-3094  (2004). 

14L.  E.  Court  and  L.  Dong,  "Automatic  registration  of  the  prostate  for 
computed-tomography-guided  radiotherapy,"  Med.  Phys.  3Q  2750-2757 
(2003). 

15E.  Schreibmann,  G.  T.  Chen,  and  L.  Xing,  "Image  interpolation  in  4D  CT 
using  a  BSpline  deformable  registration  model,"  Int.  j.  Radiat.  Oncol., 
Biol.,  Phys.  64,  1537-1550  (2006). 

16W.  Lu  etal.,  "Automatic  re-contouring  in  4D  radiotherapy,"  Phys.  Med. 
Biol.  51  1077-1099  (2006). 


Medical  Physics,  Vol.  35,  No.  10,  October  2008 


4459 


Xie  etal.\  Feature-based  rectal  contour  propagation 


4459 


17T.  Zhang  et  a/.,  “Automatic  delineation  of  on-line  head-and-neck  com¬ 
puted  tomography  images:  toward  on-line  adaptive  radiotherapy,"  Int.  J. 
Radiat.  Oncol.,  Biol.,  Phys.  68,  522-530  (2007). 

18U.  Malsch  eta/.,  "An  enhanced  block  matching  algorithm  for  fast  elastic 
registration  in  adaptive  radiotherapy,"  Phys.  Med.  Biol.  5 X  4789-4806 
(2006). 

19D.  Paquin,  L.  Xing,  and  D.  Levy,  "M  ultiscale  deformable  registration  of 
noisy  medical  images,"  Math.  Biosci.  Eng.  5,  125-144  (2008). 

20M  .  J .  M  urphy  et  a/.,  "How  does  CT  image  noise  affect  3D  deformable 
image  registration  for  image-guided  radiotherapy  planning?,"  M  ed.  Phys. 
35>  1145-1153  (2008). 

21M .  Staring,  S.  Klein,  and  J.  P.  W.  Pluim,  "A  rigidity  penalty  term  for 
nonrigid  registration,"  Med.  Phys.  34,  4098-4108  (2007). 

22E.  Schreibmann  eta/.,  "Image  interpolation  in  4D  CT  using  a  B  spline 
deformable  registration  model,"  Med.  Phys.  32,  1924  (2005). 

23R.  Alterovitz  et  a/.,  "Registration  of  M  R  prostate  images  with  biome¬ 
chanical  modeling  and  nonlinear  parameter  estimation,"  Med.  Phys.  33, 
446-454  (2006). 

24S.  Roberts  and  L.  Stals,  "Discrete  thin  plate  spline  smoothing  in  3D," 
Annal.  Chim.  4 5  C646-C659  (2004). 

25B.  Fei,  C.  Kemper,  and  D.  L.  Wilson,  "A  comparative  study  of  warping 
and  rigid  body  registration  for  the  prostate  and  pelvic  MR  volumes," 
Comput.  Med.  Imaging  Graph.  4,  267-281  (2003). 


26M  .  M  .  Coselmon  eta/.,  "M  utual  information  based  CT  registration  of  the 
lung  at  exhale  and  inhale  breathing  states  using  thin-plate  splines,"  M  ed. 
Phys.  31  2942-2948  (2004). 

27E.  Schreibmann  and  L.  Xing,  "Image  registration  with  auto-mapped  con¬ 
trol  volumes,"  Med.  Phys.  33,  1165-1179  (2006). 

28B.  Kim  eta/.,  "Mutual  information  for  automated  unwarping  of  rat  brain 
autoradiographs,"  Neuroimage  5,  31-40  (1997). 

29C.  R.  Meyer  eta/.,  "Demonstration  of  accuracy  and  clinical  versatility  of 
mutual  information  for  automatic  multimodality  image  fusion  using  affine 
and  thin-plate  spline  warped  geometric  deformations,"  Med.  I  mage  Anal. 
I  195-206  (1997). 

30K.  M  .  Brock  eta/.,  "Automated  generation  of  a  four-dimensional  model 
of  the  liver  using  warping  and  mutual  information,"  Med.  Phys.  3Q 
1128-1133  (2003). 

31L.  Xing,  J.  Siebers,  and  P.  Keall,  "Computational  challenges  for  image- 
guided  radiation  therapy:  Framework  and  current  research,"  Semin.  Ra¬ 
diat.  Oncol.  17  245-257  (2007). 

32H.  Lester  and  S.  R.  Arridge,  "A  survey  of  hierarchical  non-linear  medical 
image  registration,"  Pattern  Recogn.  32,  129-149  (1999). 

33J .  V.  Hajnal,  D.  L.  G.  Hill,  and  D.  J .  Hawkes,  Medical  Image  Registra¬ 
tion.  (CRC  Press,  Boca  Raton,  1999). 

34ITK,  http://www.itk.org. 

35VTK ,  http://public.kitware.com/VTK/. 


Medical  Physics,  Vol.  35,  No.  10,  October  2008 


IOP  Publishing 


Physics  in  Medicine  and  Biology 


Phys.  Med.  Biol.  53  (2008)  4533-4542 


doi:  10.108 8/003 1-9155/53/1 7/005 


Auto-propagation  of  contours  for  adaptive  prostate 
radiation  therapy 


Ming  Chao,  Yaoqin  Xie  and  Lei  Xing 

Department  of  Radiation  Oncology,  Stanford  University  School  of  Medicine, 

875  Blake  Wilbur  Drive,  Stanford,  CA  94305-5847,  USA 

E-mail:  lei@reyes.stanford.edu 

Received  13  March  2008,  in  final  form  9  June  2008 

Published  1  August  2008 

Online  at  stacks.iop.org/PMB/53/4533 

Abstract 

The  purpose  of  this  work  is  to  develop  an  effective  technique  to  automatically 
propagate  contours  from  planning  CT  to  cone  beam  CT  (CBCT)  to  facilitate 
CBCT-guided  prostate  adaptive  radiation  therapy.  Different  from  other  disease 
sites,  such  as  the  lungs,  the  contour  mapping  here  is  complicated  by  two  factors: 
(i)  the  physical  one-to-one  correspondence  may  not  exist  due  to  the  insertion 
or  removal  of  some  image  contents  within  the  region  of  interest  (ROI);  and  (ii) 
reduced  contrast  to  noise  ratio  of  the  CBCT  images  due  to  increased  scatter. 
To  overcome  these  issues,  we  investigate  a  strategy  of  excluding  the  regions 
with  variable  contents  by  a  careful  design  of  a  narrow  shell  signifying  the 
contour  of  an  ROI.  For  rectum,  for  example,  a  narrow  shell  with  the  delineated 
contours  as  its  interior  surface  was  constructed  to  avoid  the  adverse  influence 
of  the  day-to-day  content  change  inside  the  rectum  on  the  contour  mapping. 
The  corresponding  contours  in  the  CBCT  were  found  by  warping  the  narrow 
shell  through  the  use  of  B Spline  deformable  model.  Both  digital  phantom 
experiments  and  clinical  case  testing  were  carried  out  to  validate  the  proposed 
ROI  mapping  method.  It  was  found  that  the  approach  was  able  to  reliably  warp 
the  constructed  narrow  band  with  an  accuracy  better  than  1.3  mm.  For  all  five 
clinical  cases  enrolled  in  this  study,  the  method  yielded  satisfactory  results  even 
when  there  were  significant  rectal  content  changes  between  the  planning  CT 
and  CBCT  scans.  The  overlapped  area  of  the  auto-mapped  contours  over  90% 
to  the  manually  drawn  contours  is  readily  achievable.  The  proposed  approach 
permits  us  to  take  advantage  of  the  regional  calculation  algorithm  yet  avoiding 
the  nuisance  of  rectum/bladder  filling  and  provide  a  useful  tool  for  adaptive 
radiotherapy  of  prostate  in  the  future. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 
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1.  Introduction 

There  are  increased  interests  in  online/off-line  adaptive  replanning  with  consideration  of 
updated  patient  anatomy  information  derived  from  cone  beam  CT  (CBCT)  (Yan  et  al  2000, 
Court  et  al  2005,  Ghilezan  et  al  2005,  Langen  et  al  2005,  Letourneau  et  al  2005,  Oldham 
et  al  2005,  Xing  et  al  2006,  Yang  et  al  2007).  An  indispensable  step  toward  the  realization 
of  this  type  of  adaptive  replanning  is  the  segmentation  of  CBCT  images  (Gao  et  al  2006, 
Smeenk  et  al  2007).  While  this  task  is,  in  principle,  achievable  using  deformable  registration 
of  the  planning  CT  and  CBCT  images,  the  accuracy  of  the  registration,  and  therefore  contour 
mapping,  is  often  adversely  affected  by  the  presence  of  image  contents  in  one  image  that  do  not 
have  correspondence  in  the  other  image,  such  as  the  presence/absence  of  bowel  gas  or  bladder 
filling  in  the  pelvic  patients  (Huang  et  al  2002).  A  regional  contour  mapping  algorithm  is  highly 
desirable  for  its  high  computational  efficiency  and  accuracy  (Schreibmann  and  Xing  2005, 
Chao  et  al  2007,  2008).  It  also  has  the  potential  to  improve  the  robustness  of  deformable 
registration  because  the  mapped  contours  can  be  utilized  as  a  priori  system  information  to 
facilitate  the  registration  of  the  two  input  images. 

Large  interfraction  organ  motion  exists  in  radiation  therapy  of  pelvic  diseases,  such  as 
prostate  and  cervical  cancers  (Balter  et  al  1995,  Yan  and  Lockman  2001,  Michalski  2003). 
Adaptive  therapy  provides  a  viable  option  to  ensure  an  adequate  dose  coverage  of  the  tumor 
target  while  sparing  the  sensitive  structures  (Yan  et  al  1997,  Court  et  al  2005,  Wu  et  al  2008). 
Toward  the  general  goal  of  clinical  implementation  of  adaptive  therapy  of  prostate  cancer, 
in  this  work  we  develop  an  effective  regional  contour  mapping  technique  to  automatically 
propagate  rectum  and  prostate  contours  from  planning  CT  to  CBCT.  Different  from  other 
organs,  such  as  the  lungs  and  liver,  the  contour  mapping  here  is  complicated  by  two  factors: 

(1)  the  physical  one-to-one  correspondence  may  not  exist  due  to  the  insertion  or  removal  of 
some  image  contents  within  the  rectum  (Gao  et  al  2006);  and  (ii)  reduced  contrast  to  noise 
ratio  (CNR)  of  the  CBCT  images  due  to  increased  scatter  in  CBCT  scanning  (Li  et  al  2007). 
A  solution  customized  to  rectum  mapping  is  proposed  to  overcome  the  above  two  limitations 
and  allow  us  to  take  advantage  of  the  regional  calculation  algorithm  yet  avoiding  the  nuisance 
of  rectum/bladder  filling. 

2.  Methods  and  materials 

2.7.  Software  platform 

The  narrow  band  contour  mapping  algorithm  was  implemented  using  the  NLM  Insight  Toolkit 
(ITK)  (Ibanez  et  al  2003)  and  the  Visualization  Toolkit  (VTK)  (Schroeder  et  al),  which  are 
open  source  cross-platform  C++  software  toolkits.  They  are  freely  available  for  research 
purposes  (http://www.itk.org  for  ITK  and  http://public.kitware.com/VTK/  for  VTK).  ITK 
provides  various  kinds  of  basic  algorithms  for  performing  registration  and  segmentation  for 
medical  images.  The  fact  that  programs  contained  in  ITK  are  highly  extendable  makes  it  an 
ideal  platform  for  the  development  of  image  registration  methods.  VTK  is  primarily  used  for 
image  visualization. 

2.2.  The  contour  propagation  process 

The  rectum  contour  propagation  from  planning  CT  to  CBCT  was  carried  out  in  the  following 
sequences:  (1)  rigid  alignment  of  the  planning  CT  and  CBCT  images  using  a  rigid  registration; 

(2)  construction  of  a  narrow  band  extended  to  a  region  1-2  cm  outside  the  manually  delineated 
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Narrow  Shell 
on  Planning  CT 


Narrow  Shell  Warping 
on  CBCT 


Mapped  Contour  on  CBCT 


Figure  1.  Flowchart  of  the  narrow  shell  contour  mapping  process. 


Figure  2.  Axial  view  (a),  coronal  view  (b),  sagittal  view  (c)  and  volumetric  rendering  (d)  of 
planning  CT  of  the  first  patient.  Template  contours  of  prostate  (magenta)  and  rectum  (green)  were 
superimposed  on  CT  images  as  indicated. 


contours  on  planning  CT;  (3)  warping  of  the  shell  from  the  planning  CT  to  CBCT.  Upon 
successful  warping  of  the  shell,  the  deformation  field  in  step  (3)  was  utilized  to  transform  the 
manually  segmented  contours  to  the  CBCT  images.  The  central  idea  here  is  to  exclude  the 
volume  inside  the  rectum  because  its  content  may  change  from  day  to  day.  Prostate  contour 
can  be  mapped  similarly,  but  the  shell  spans  the  regions  both  inside  and  outside  the  manually 
segmented  prostate.  Figure  1  depicts  the  flowchart  of  the  overall  contour  mapping  process 
using  the  narrow  shell  technique. 

The  manually  delineated  region  of  interest  (ROI)  contours  on  planning  CT  are  referred  to 
as  the  template  contours.  Figure  2  illustrates  the  planning  CT  and  the  template  contours  for 
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Figure  3.  Narrow  shell  representation  of  regions  of  interest  (ROIs).  The  curves  represent  contours 
and  the  shaded  regions  sketch  the  narrow  shell  construct.  Illustrated  are  axial  view  (a),  coronal  view 
(b)  and  sagittal  view  (d)  of  the  planning  CT  image  with  a  single-sided  narrow  shell  surrounding 
the  rectum  template  contour. 


the  ROIs  such  as  prostate  and  rectum  for  one  of  the  patients  studied  in  this  work.  A  narrow 
shell  was  constructed  using  the  method  described  above  for  each  structure.  Typical  shapes 
of  the  narrow  shells  for  rectum  are  illustrated  in  figure  3.  The  shell  width  was  chosen  to  be 
10  mm  in  this  study. 

Upon  the  construction  of  shell  structure,  the  shell  is  warped  onto  CBCT  images  by  using 
a  B Spline  model  (Lee  et  al  1996,  1997,  Badea  et  al  2008)  with  a  normalized  cross  correlation 
metric  function  defined  as 

y/lliLi  (J/V))2' 

where  /A(x)  and  IB(x') are  the  intensities  of  the  ith  voxels  of  images  A  and  B,  respectively; 
x'  =  Tx,  where  T  is  the  transformation  matrix.  In  the  narrow  shell  warping  calculation,  the 
input  images  to  the  registration  software  are  the  narrow  shell,  /A(x),  and  the  CBCT  image, 
IB(x').  The  narrow  shell  acts  as  a  shape  representation  model  of  an  anatomic  structure.  The 
task  here  is  to  find  the  deformation,  T(x),  that  maps  an  arbitrary  point  v  on  the  fixed  image  to 
the  corresponding  point  x  on  the  moving  image  (or  vice  versa)  so  that  the  best  possible  match, 
as  measured  by  the  registration  metric,  is  achieved.  The  calculation  proceeds  in  an  iterative 
fashion.  The  limited  memory  Broyden-Fletcher-Goldfarb-Shannon  algorithm  (L-BFGS)  (Liu 
and  Nocedal  1989)  was  used  for  the  iterative  optimization.  The  L-BFGS  optimizer  is  well 
known  for  its  superior  performance  in  dealing  with  high  dimensional  problems  (Schreibmann 
and  Xing  2005,  Schreibmann  et  al  2006).  We  briefly  describe  the  optimization  algorithm  in 
the  following. 

The  BFGS  method  is  derived  from  Newton’s  method  in  optimization  which  is  a  class  of 
hill-climbing  optimization  techniques.  It  tries  to  seek  the  stationary  point  of  a  function,  where 
the  gradient  is  0.  Starting  from  a  positive  definitive  approximation  of  the  inverse  Hessian  H0 
at  x0,  L-BFGS  derives  the  optimization  variables  by  iteratively  searching  through  the  solution 
space.  At  an  iteration  k,  the  calculation  proceeds  as  follows: 

(1)  determine  the  descent  direction  =  —  HkS7f(xk); 

(2)  line  search  with  a  step  size  otk  =  arga^0  min  f(xk  +  ap&),  where  a  is  the  step  size  defined 
in  the  L-BFGS  software  package; 

(3)  update  x^+i  =  xk  +  ak  p^; 

(4)  compute  H^+i  with  the  updated  H^. 

At  each  iteration  a  backtracking  line  search  is  used  in  L-BFGS  to  determine  the  step  size  of 
movement  to  reach  the  minimum  of/along  the  ray  xk  +  apk.  During  optimization,  the  iterative 
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(a)  (b) 


Figure  4.  Digital  phantom  experiment,  (a)  Original  CT  image  and  manual  contour;  the  vectors 
overlaid  on  the  image  are  the  displacement  field  introduced  to  deform  the  image,  (b)  Digitally 
deformed  CT  image  with  ground  truth  contours  and  mapped  contours.  The  template  contour  after 
rigid  registration  of  the  two  images  is  also  displayed. 


calculation  continues  until  the  following  stopping  criterion  is  fulfilled: 

I|V/(X*)||2 

-  <  £, 

max(l,  ||x*||2) 

or  a  pre-set  maximum  number  of  iterations  is  reached.  In  this  study  we  set  e  =  10-6  and  the 
iteration  number  to  200. 

2.3.  Evaluation  of  algorithm  performance 

Evaluation  of  a  contour  propagation  algorithm  is  often  a  difficult  task  because  of  the  lack  of  the 
ground  truth  for  comparison.  A  commonly  used  approach  is  to  visually  inspect  the  mapped 
contours  by  a  clinician.  The  use  of  synthetic  images  (digital  phantoms)  is  also  useful.  In 
this  test,  the  images  and  existing  contours  are  distorted  with  known  deformation  fields,  which 
serves  as  the  ground  truth  for  quantitative  evaluation  of  the  model-based  contour  propagation. 
In  constructing  the  digital  phantom,  the  deformation  field  used  to  deform  a  planning  pelvic 
CT  image  was  generated  as  follows.  B Spline  8x8x8  grids  were  overlaid  on  the  CT  image. 
Three  random  generators  following  the  Gaussian  distribution  were  called  for  to  provide  the 
displacements  for  the  components  of  x,  y  and  z  at  each  grid  point.  The  mean  value  /x  of  the 
Gaussian  function  were  all  chosen  to  be  0,  but  the  variances  of  the  three  components  were  set 
to  be 

=  cfy  =  5.0  mm 
<7^  ~  2.0  mm. 

Figure  4(a)  shows  the  deformation  field  overlaid  on  the  original  CT  image.  The  manual 
contour  on  the  original  CT  image  was  warped  with  the  same  deformation  field.  This  contour 
served  as  the  ground  truth.  A  comparison  of  the  mapped  contour  with  the  ground  truth  allowed 
us  to  quantitatively  assess  the  accuracy  of  the  proposed  approach. 

2.4.  Patient  study 

The  image  data  sets  for  five  prostate  cancer  patients  were  used  to  further  evaluate  the  proposed 
technique.  The  CBCT  images  were  acquired  using  a  Trilogy™  system  (Varian  Medical 
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Figure  5.  Rectal  contours  on  axial,  coronal  and  sagittal  planes  obtained  with  the  proposed  narrow 
shell-based  (in  red)  and  conventional  B Spline  approach  (in  yellow)  for  the  first  clinical  case.  The 
contours  after  rigid  alignment  (in  blue)  were  overlaid  on  the  same  CT  images. 


Systems,  Palo  Alto,  CA).  All  patient  CBCT  images  were  acquired  with  90  mA  tube  current 
and  120  kV  voltage.  Under  this  image  acquisition  protocol,  the  dose  was  around  4cGy  for 
a  kV  CBCT  scan  (Wen  et  al  2007).  The  CBCT  image  sets  for  all  enrolled  patients  were 
reconstructed  with  a  2.5  mm  slice  thickness.  The  planning  CT  scans  were  acquired  with  a  GE 
Discovery-ST  CT  scanner  (GE  Medical  System,  Milwaukee,  WI).  The  size  of  a  CT  slice  was 
512  x  512.  The  planning  CT  images  were  reconstructed  with  a  1.25  mm  slice  thickness  for 
all  patients.  All  CT  images  were  transferred  through  DICOM  to  a  high-performance  personal 
computer  (PC)  with  Xeon  (3.6  GHz)  processor  for  image  processing. 


3.  Results  and  discussions 

3.1.  Digital  phantom  study 

The  narrow  shell  constructed  based  on  the  planning  CT  image  and  the  template  rectal  contour 
(figure  4(a))  was  warped  to  the  digitally  deformed  images  using  the  technique  described  above. 
For  convenience,  the  deformation  field  used  to  deform  the  planning  CT  images  is  displayed  in 
figure  4(a)  as  a  vector  map.  The  rectal  contour  was  then  extracted  upon  the  successful  warping 
of  the  shell.  The  resultant  rectum  contour  is  shown  in  figure  4(b)  together  with  the  ground 
truth  contour.  The  difference  between  the  two  sets  of  contours  was  found  to  be  insignificant 
by  visual  inspection.  A  quantitative  comparison  indicated  that  the  mean  discrepancy  between 
the  two  sets  of  contours  was  less  than  1.3  mm. 

3.2.  Clinical  case  study 

The  proposed  approach  was  applied  to  five  enrolled  prostate  cancer  patients.  Figure  5  shows 
the  result  of  contour  propagation  for  the  first  clinical  case  in  axial,  coronal  and  sagittal  views 
of  CBCT  image,  respectively.  For  comparison,  the  result  obtained  by  registering  the  whole 
images  using  the  conventional  B  Spline  method  was  also  superimposed  in  figure  5.  We  also 
observed  the  convergence  behaviors  of  the  two  approaches  by  examining  the  metric  function 
as  a  function  of  the  iteration  step.  Nearly  30  iterations  were  needed  in  the  conventional 
registration,  whereas  only  12  were  necessary  to  obtain  converged  result  in  the  narrow  shell- 
based  approach.  In  this  particular  case,  both  approaches  yielded  acceptable  results  because 
of  no  substantial  change  in  rectal  filling  in  planning  CT  and  CBCT.  The  difference  between 
two  sets  of  mapping  is  well  within  2  mm.  However,  the  narrow  shell-based  calculation  was 
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Figure  6.  Propagated  rectal  contour  (in  red)  overlaid  on  CBCT  images  in  the  second  patient  case. 
The  results  based  on  conventional  approach  (in  yellow)  and  physician’s  manual  contour  were  also 
superimposed  on  CBCT  images  as  indicated. 


Figure  7.  Contour  propagation  results  for  the  fourth  clinical  case,  (a),  (b)  and  (c)  The  rectum 
contours  on  axial,  coronal  and  sagittal  views  of  CBCT  image;  (d),  (e)  and  (f)  the  mapped  prostate 
contours  on  CBCT  images.  The  manual  contours  were  also  overlaid  on  the  same  CBCT  images. 


found  to  be  an  order  of  magnitude  more  efficient  as  measured  by  the  computational  speed,  in 
addition  to  the  greatly  reduced  usage  of  computer  memory.  Results  for  the  second  patient  are 
shown  in  figure  6.  For  this  patient,  however,  the  rectal  fillings  in  planning  CT  and  CBCT  were 
different  and  the  conventional  B Spline  approach  led  to  a  larger  error  (mean  error  ~8.0  mm) 
or  failed  to  accurately  map  the  rectum  from  planning  CT  to  CBCT.  The  physician’s  manual 
contour  was  also  superimposed  on  the  images  for  visual  comparison. 

The  mapped  rectum  and  prostate  contours  for  another  case  are  displayed  in  figure  7. 
Overall,  the  algorithm  performed  well  for  both  rectum  and  prostate.  The  mean  error  of  the 
mapped  contours  was  estimated  to  be  around  2  mm  for  the  rectum,  and  2.5  mm  for  the  prostate. 
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The  slightly  larger  error  in  prostate  contour  mapping  is  attributed  to  the  fact  that  the  image 
contrast  in  the  prostate  boundary  region  was  not  high.  In  this  case,  a  careful  examination  and 
even  a  manual  modification  of  the  mapped  contour  by  a  physician  may  be  needed  in  clinical 
operation. 

Up  to  this  point,  little  has  been  done  to  register  images  with  the  general  problem  of  image 
content  change  as  studied  in  this  work.  Foskey  et  al  (2005)  and  Gao  et  al  (2006)  studied  the 
topics  by  using  different  models.  The  essence  of  their  work  was  to  introduce  ‘artificial  gas’  to 
mimic  the  rectal  image  content  or  shrink  the  ‘gas’  region  by  deflation.  While  the  introduction 
of  the  artificial  change  could  improve  the  registration  accuracy  in  some  special  cases,  the 
approaches  may  fail  when  the  change  in  rectum  content  involves  solid  substances,  such  as  the 
bowel.  A  narrow  shell  model  seems  to  be  natural  in  addressing  the  issue  of  image  content 
inconsistence.  The  technique  is  conceptually  simple,  easy  to  implement  and  valid  for  a  wide 
variety  of  clinical  situations. 

One  of  the  practical  concerns  is  that  the  relatively  low  quality  of  CBCT  images  may 
influence  the  accuracy  of  image  registration  and  thus  the  contour  mapping.  Paquin  et  al  (2008) 
quantitatively  studied  the  influence  of  different  types  of  noises  on  deformable  registration  and 
found  that  the  accuracy  of  image  registration  does  not  depend  on  the  global  noise  unless  the 
noise  reaches  a  certain  threshold  value.  Murphy  et  al  (2008)  also  demonstrated  that  noise 
levels  in  cone-beam  CTs  that  might  reduce  manual  contouring  accuracy  do  not  reduce  image 
registration  and  automatic  contouring  accuracy. 

In  principal,  the  proposed  method  is  also  applicable  to  facilitate  the  contour  warping  from 
planning  CT  to  MV  CBCT  or  images  acquired  using  the  tomotherapy  devices.  However,  it 
is  important  to  emphasize  that  a  prerequisite  of  any  deformable  model-based  approach  is  that 
there  are  sufficient  image  features  in  the  neighborhood  of  the  contours  to  guide  the  warping  of 
the  images.  Because  of  the  inherent  low  contrast  of  the  MV  images,  the  deformable  registration 
of  CT  and  MV  CBCT  may  be  challenging.  Lu  et  al  (2004)  have  reported  successful  deformable 
registration  of  CT  and  MV  CT  images.  It  will  be  an  interesting  subject  of  future  research 
to  examine  whether  a  similar  or  better  accuracy  of  registration  can  be  achieved  when  a  local 
deformable  registration  technique  is  used  for  the  system. 


4.  Summary 

An  automatic  contour  mapping  technique  customized  for  adaptive  radiation  therapy  of  prostate 
cancer  has  been  described.  Deformable  image  registration  has  been  extensively  studied.  The 
narrow  shell  constructed  surrounding  a  ROI  contour  provides  a  compact  representation  of 
the  ROI  and  greatly  facilitates  the  contour  mapping  calculation.  The  chief  advantage  of  the 
proposed  technique  is  that  it  is  much  less  prone  to  the  variation  in  the  image  contents,  which 
occurs  frequently  in  the  pelvic  region  due  to  involuntary  physiological  process.  Additionally, 
the  approach  is  computationally  more  efficient  with  lower  computer  memory  consumption  and 
fast  computational  speed.  Both  phantom  and  clinical  case  studies  yielded  satisfactory  results 
and  suggested  that  the  method  may  prove  to  be  useful  for  adaptive  radiotherapy  of  prostate 
cancer  in  the  future. 
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